﻿using AEFunction;
using AETools.UI;
using CommonFunction;
using DbUtility;
using CloneClass;
using ESRI.ArcGIS.Carto;
using ESRI.ArcGIS.Catalog;
using ESRI.ArcGIS.CatalogUI;
using ESRI.ArcGIS.Controls;
using ESRI.ArcGIS.DataSourcesRaster;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.SystemUI;
using System;
using System.Collections.Generic;
using System.Data;
using System.Data.OleDb;
using System.Diagnostics;
using System.IO;
using System.Windows.Forms;
using System.Xml;
using ESRI.ArcGIS.Geoprocessor;
using System.Collections;
using System.Text.RegularExpressions;
using ESRI.ArcGIS.DataSourcesFile;
using ESRI.ArcGIS.DataSourcesGDB;
using AEFunction.UI;
namespace AETools
{

    public sealed partial class MainForm : Form
    {
        #region class private members
        private IMapControl3 m_mapControl = null;
        private string m_mapDocumentName = string.Empty;
        private AEFunction.UI.AttributeTableDlg attributeDlg = new AEFunction.UI.AttributeTableDlg();
        private string rulePath = Application.StartupPath + "\\rules.mdb";
        private AccessOperation ao = null;

        #region 元胞自动机变量

        IRasterLayer pFirstRasterLayer = null;//第一次影像
        IRasterLayer pSecRasterLayer = null;//第二次影像
        IRasterLayer pForecastRasterLayer = null;//模拟影像
        Dictionary<int, LandType> ltList = new Dictionary<int, LandType>();//地类编码与地类类型的对应表
        Dictionary<int, int> dlbmToIndex = new Dictionary<int, int>();//地类编码的顺序
        Dictionary<int, int> IndexTodlbm = new Dictionary<int, int>();//根据顺序得到地类编码
        IFeatureLayer pForecastFeatureLayer = null;//模拟影像对应的点图层
        //Dictionary<LandType, List<Land>> allLand = new Dictionary<LandType, List<Land>>();//按地类类型分类所有地块
        //Dictionary<LandType, List<Land>> forecastAllLand = new Dictionary<LandType, List<Land>>();
        Dictionary<int, Dictionary<int, LogisticFactor>> LogisticsFactors = new Dictionary<int, Dictionary<int, LogisticFactor>>();//求得的逻辑回归参数
        List<int> fjsyd = new List<int>();//非建设用地
        //Dictionary<Position, Land> posLand = new Dictionary<Position, Land>();//按行和列记录地块，方便在元胞模拟是处理
        int[] cellValue = null;
        string directoryPath = "";
        #endregion
        #endregion

        #region class constructor
        public MainForm()
        {
            InitializeComponent();
        }
        #endregion

        private void MainForm_Load(object sender, EventArgs e)
        {
            //get the MapControl
            m_mapControl = (IMapControl3)axMapControl1.Object;

            //disable the Save menu (since there is no document yet)
            menuSaveDoc.Enabled = false;
            ao = new AccessOperation(rulePath);
            fjsyd.Add(1);
            fjsyd.Add(2);
            fjsyd.Add(3);
            fjsyd.Add(4);
            fjsyd.Add(6);
            fjsyd.Add(7);
            //fjsyd.Add(8);
        }

        #region Main Menu event handlers
        private void menuNewDoc_Click(object sender, EventArgs e)
        {
            //execute New Document command
            ICommand command = new CreateNewDocument();
            command.OnCreate(m_mapControl.Object);
            command.OnClick();
        }

        private void menuOpenDoc_Click(object sender, EventArgs e)
        {
            //execute Open Document command
            ICommand command = new ControlsOpenDocCommandClass();
            command.OnCreate(m_mapControl.Object);
            command.OnClick();
        }

        private void menuSaveDoc_Click(object sender, EventArgs e)
        {
            //execute Save Document command
            if (m_mapControl.CheckMxFile(m_mapDocumentName))
            {
                //create a new instance of a MapDocument
                IMapDocument mapDoc = new MapDocumentClass();
                mapDoc.Open(m_mapDocumentName, string.Empty);

                //Make sure that the MapDocument is not readonly
                if (mapDoc.get_IsReadOnly(m_mapDocumentName))
                {
                    MessageBox.Show("Map document is read only!");
                    mapDoc.Close();
                    return;
                }

                //Replace its contents with the current map
                mapDoc.ReplaceContents((IMxdContents)m_mapControl.Map);

                //save the MapDocument in order to persist it
                mapDoc.Save(mapDoc.UsesRelativePaths, false);

                //close the MapDocument
                mapDoc.Close();
            }
        }

        private void menuSaveAs_Click(object sender, EventArgs e)
        {
            //execute SaveAs Document command
            ICommand command = new ControlsSaveAsDocCommandClass();
            command.OnCreate(m_mapControl.Object);
            command.OnClick();
        }

        private void menuExitApp_Click(object sender, EventArgs e)
        {
            //exit the application
            Application.Exit();
        }
        #endregion

        //listen to MapReplaced evant in order to update the statusbar and the Save menu
        private void axMapControl1_OnMapReplaced(object sender, IMapControlEvents2_OnMapReplacedEvent e)
        {
            //get the current document name from the MapControl
            m_mapDocumentName = m_mapControl.DocumentFilename;

            //if there is no MapDocument, diable the Save menu and clear the statusbar
            if (m_mapDocumentName == string.Empty)
            {
                menuSaveDoc.Enabled = false;
                statusBarXY.Text = string.Empty;
            }
            else
            {
                //enable the Save manu and write the doc name to the statusbar
                menuSaveDoc.Enabled = true;
                statusBarXY.Text = System.IO.Path.GetFileName(m_mapDocumentName);
            }
        }

        private void axMapControl1_OnMouseMove(object sender, IMapControlEvents2_OnMouseMoveEvent e)
        {
            statusBarXY.Text = string.Format("{0}, {1}  {2}", e.mapX.ToString("#######.##"), e.mapY.ToString("#######.##"), axMapControl1.MapUnits.ToString().Substring(4));
        }

        private void readExcelGetFeatureToolStripMenuItem_Click(object sender, EventArgs e)
        {
            OpenFileDialog ofd = new OpenFileDialog();
            ofd.Filter = "Excel2007(*.xlsx)|*.xlsx|Excel2003(*.xls)|*.xls";
            if (ofd.ShowDialog() == DialogResult.OK)
            {
                object missing = System.Reflection.Missing.Value;

                string excelFileName = ofd.FileName;
                string strCon = " Provider =Microsoft.ACE.OleDb.12.0 ; Data Source =" + excelFileName + ";Extended Properties='Excel 12.0; HDR=No; IMEX=1'";
                OleDbConnection myConn = new OleDbConnection(strCon);
                string strCom = "select * from [Sheet1$]";
                OleDbDataAdapter myCommand = new OleDbDataAdapter(strCom, myConn);
                myConn.Open();
                DataTable myDataTable = new DataTable();
                myCommand.Fill(myDataTable);
                ChooseMatchFields CMF1 = new ChooseMatchFields(ArcEngineFunction.GetAllLayer(axMapControl1), Function.GetExcelField(myDataTable));
                if (CMF1.ShowDialog() != DialogResult.OK)
                    return;
                //shape的匹配字段和Excel的匹配字段
                string shapeField = CMF1.ShapeField;
                string excelField = CMF1.ExcelField;
                //匹配的图层
                string layerName = CMF1.layerName;
                ILayer selectLayer = ArcEngineFunction.GetLayerByName(axMapControl1, layerName);//筛选图层

                if (selectLayer == null)
                    return;

                SelectNewShapeFields SNSF = new SelectNewShapeFields(ArcEngineFunction.GetAllLayer(axMapControl1), Function.GetExcelField(myDataTable));
                if (SNSF.ShowDialog() != DialogResult.OK)
                    return;
                //需要添加的字段
                string[] ShapeSelectedFields = SNSF.ShapeSelectedFields;
                List<ExcelField> ExcelSelectedFields = SNSF.ExcelSelectedFields;
                IFields pFields = ArcEngineFunction.GetFields(selectLayer, ShapeSelectedFields);

                string savePath = SNSF.SavePath;
                //获取文件名，不带路径
                string fileName = savePath.Substring(savePath.LastIndexOf("\\") + 1, savePath.LastIndexOf(".") - savePath.LastIndexOf("\\") - 1);
                //获取文件路径，不带文件名
                string FilePath = savePath.Substring(0, savePath.LastIndexOf("\\"));
                ChangeFieldType CFT = new ChangeFieldType(ExcelSelectedFields);
                if (CFT.ShowDialog() != DialogResult.OK)
                    return;
                ExcelSelectedFields = CFT.typeList;
                ArcEngineFunction.CreatePointShapefile(ref fileName, FilePath, pFields, ExcelSelectedFields);

                axMapControl1.AddShapeFile(FilePath, fileName + ".shp");
                ILayer newLayer = axMapControl1.get_Layer(0);
                //选择的图层字段在老图层对应的列号
                List<int> LayerFieldsFromOldLayer = ArcEngineFunction.GetFieldIndexFromLayer(selectLayer, ShapeSelectedFields);
                //选择的图层字段在新图层对应的列号
                List<int> LayerFieldsFromNewLayer = ArcEngineFunction.GetFieldIndexFromLayer(newLayer, ShapeSelectedFields);
                //选择的Excel字段在Excel对应的列号
                List<int> ExcelFieldsFromExcel = Function.GetFieldIndexFromTable(myDataTable, ExcelSelectedFields);
                //选择的Excel字段在新图层对应的列号
                List<int> ExcelFieldsFromNewLayer = ArcEngineFunction.GetExcelFieldIndexFromLayer(newLayer, ExcelSelectedFields.ToArray());

                IFeatureClass fc = AEFunction.ArcEngineFunction.GetLayerFeatureClass(newLayer);
                IDataset dataset = (IDataset)fc;
                IWorkspace workspace = dataset.Workspace;
                //Cast for an IWorkspaceEdit  
                IWorkspaceEdit workspaceEdit = (IWorkspaceEdit)workspace;
                //start an edit session and operation              
                workspaceEdit.StartEditing(true);
                workspaceEdit.StartEditOperation();
                int index = Function.GetFieldIndex(myDataTable, excelField);
                string nullstr = "";
                for (int i = 1; i < myDataTable.Rows.Count; i++)
                {
                    #region 构建选择语句
                    string str = myDataTable.Rows[i][index].ToString().Trim();
                    IQueryFilter queryfilter = new QueryFilterClass();
                    queryfilter.SubFields = shapeField;
                    queryfilter.WhereClause = String.Format("{0}='{1}'", shapeField, str);
                    #endregion
                    IFeatureCursor featureCursor = ArcEngineFunction.GetLayerFeatureClass(selectLayer).Search(queryfilter, true);
                    IFeature feature = featureCursor.NextFeature();
                    if (feature == null)
                    {
                        nullstr += str + ",";
                    }
                    while (feature != null)
                    {
                        IPoint p = ArcEngineFunction.GetLayerFeatureClass(selectLayer).GetFeature(feature.OID).Shape as IPoint;
                        IFeature f = fc.CreateFeature();
                        f.Shape = p as IGeometry;
                        #region 属性赋值
                        for (int j = 0; j < LayerFieldsFromOldLayer.Count; j++)
                        {
                            f.set_Value(LayerFieldsFromNewLayer[j], feature.get_Value(LayerFieldsFromOldLayer[j]));
                        }

                        for (int j = 0; j < ExcelFieldsFromExcel.Count; j++)
                        {
                            string tmp = myDataTable.Rows[i][ExcelFieldsFromExcel[j]].ToString().Trim();
                            if (ExcelSelectedFields[j].type == Type.GetType("System.Int32") || ExcelSelectedFields[j].type == Type.GetType("System.Double"))
                                if (tmp == "")
                                    tmp = "0";
                            f.set_Value(ExcelFieldsFromNewLayer[j], (object)tmp);
                        }
                        #endregion
                        f.Store();
                        feature = featureCursor.NextFeature();
                    }
                }

                workspaceEdit.StopEditOperation();
                workspaceEdit.StopEditing(true);
                IFeatureClassManage featureClassManage = (IFeatureClassManage)fc;
                featureClassManage.UpdateExtent();
                MessageBox.Show("提取完毕。", "提示", MessageBoxButtons.OK);
            }
        }

        private void axTOCControl1_OnMouseDown(object sender, ITOCControlEvents_OnMouseDownEvent e)
        {
            if (axMapControl1.LayerCount <= 0)
                return;
            esriTOCControlItem pTOCControlItem = new esriTOCControlItem();
            IBasicMap pBasicMap = new MapClass();
            ILayer pLayer = new FeatureLayerClass();
            object pOther = new object();
            object pIndex = new object();
            axTOCControl1.HitTest(e.x, e.y, ref pTOCControlItem,
                ref pBasicMap, ref pLayer, ref pOther, ref pIndex);
            if (e.button == 2)
            {
                if (pBasicMap == null)
                    return;
                CMenuStrip cms = new CMenuStrip(pLayer, axMapControl1, ref attributeDlg);
                cms.ContextMenuStrip.Show(axTOCControl1, e.x, e.y);
            }
        }

        private void axMapControl1_OnMouseDown(object sender, IMapControlEvents2_OnMouseDownEvent e)
        {
            switch (e.button)
            {
                case 4://中键
                    {
                        this.Cursor = Cursors.NoMove2D;
                        axMapControl1.Pan();

                        break;
                    }
                case 2://中键
                    {
                        axMapControl1.CurrentTool = null;
                        break;
                    }
            }
        }

        private void ExportScreenXY_Click(object sender, EventArgs e)
        {
            List<ILayer> allLayer = new List<ILayer>();
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                allLayer.Add(axMapControl1.get_Layer(i));
            }
            SelectLayers sl = new SelectLayers(allLayer.ToArray());
            if (DialogResult.OK == sl.ShowDialog())
            {
                XmlDocument xmlDoc = new XmlDocument();
                XmlElement xmlOriginalData = xmlDoc.CreateElement("OriginalData");
                xmlDoc.AppendChild(xmlOriginalData);
                foreach (ILayer l in sl.Layers)
                {
                    ExportLayer(l, xmlDoc, xmlOriginalData, m_mapControl.ActiveView);
                }
                xmlDoc.Save(sl.SavePath);
                MessageBox.Show("保存成功！");
            }
        }

        private void ExportLayer(ILayer layer, XmlDocument xmlDoc, XmlElement xmlElement, IActiveView pActiveView)
        {
            IFeatureLayer pLayer = (IFeatureLayer)layer;
            IFeatureCursor pCursor;
            pCursor = pLayer.Search(null, false);
            IFeature pFeature;
            pFeature = pCursor.NextFeature();
            while (pFeature != null)
            {
                ExportFeature(pFeature, xmlDoc, xmlElement, pActiveView);
                pFeature = pCursor.NextFeature();
            }

        }

        private void ExportFeature(IFeature pFeature, XmlDocument xmlDoc, XmlElement xmlElement, IActiveView pActiveView)
        {
            int index = pFeature.Table.FindField("GB");
            string GBCode = Convert.ToString(pFeature.get_Value(index));
            string sqlStr = "Select SymbolCode from FeatureSymbol where FeatureCode=\"" + GBCode + "\"";
            DataTable dt = ao.ExecuteDataTable(sqlStr);
            if (dt.Rows.Count == 0)
                return;
            string symbolCode = Convert.ToString(dt.Rows[0][0]);

            XmlElement xmlFeature = xmlDoc.CreateElement("Feature");
            xmlElement.AppendChild(xmlFeature);
            ExportGeometry(pFeature, xmlDoc, xmlFeature, pActiveView);
            XmlElement xmlType = xmlDoc.CreateElement("Type");


            xmlFeature.AppendChild(xmlType);
            if (pFeature.Shape.GeometryType == esriGeometryType.esriGeometryPoint)
            {
                xmlType.InnerText = "Point";
            }
            else if (pFeature.Shape.GeometryType == esriGeometryType.esriGeometryPolyline)
            {
                xmlType.InnerText = "Line";
            }
            else
                xmlType.InnerText = "Polygon";

            XmlElement xmlCode = xmlDoc.CreateElement("SymbolCode");
            xmlFeature.AppendChild(xmlCode);
            xmlCode.InnerText = symbolCode;
        }

        private void ExportGeometry(IFeature pFeature, XmlDocument xmlDoc, XmlElement xmlElement, IActiveView pActiveView)
        {
            XmlElement xmlGeometry = xmlDoc.CreateElement("Geometry");
            xmlElement.AppendChild(xmlGeometry);
            string geo = "";
            int x = 0; int y = 0;
            if (pFeature.Shape.GeometryType.ToString() == "esriGeometryPoint")
            {
                IPoint pPoint = pFeature.Shape as IPoint;

                pActiveView.ScreenDisplay.DisplayTransformation.FromMapPoint(pPoint, out x, out y);
                geo += "(" + x + "," + y + ");";
            }
            else
            {
                IPointCollection pPointCollection = pFeature.Shape as IPointCollection;
                for (int j = 0; j < pPointCollection.PointCount; j++)
                {
                    IPoint pPoint = pPointCollection.get_Point(j);
                    pActiveView.ScreenDisplay.DisplayTransformation.FromMapPoint(pPoint, out x, out y);
                    geo += "(" + x + "," + y + ");";
                }
            }
            xmlGeometry.InnerText = geo;
        }

        private void 导出矢量瓦片ToolStripMenuItem_Click(object sender, EventArgs e)
        {

            AEFunction.UI.frmCreateTiles4MultiScale exportTiles = new AEFunction.UI.frmCreateTiles4MultiScale(this.axMapControl1);
            exportTiles.ShowDialog();
        }

        private void batchClipToolStripMenuItem_Click(object sender, EventArgs e)
        {
            AEFunction.UI.BatchClip batchClip = new AEFunction.UI.BatchClip(toolStripProgressBar1, axMapControl1);
            batchClip.ShowDialog();
        }

        private void 加载数据ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            IGxDialog openRaterDLg = new GxDialogClass();
            openRaterDLg.AllowMultiSelect = true;
            openRaterDLg.Title = "打开文件";

            IEnumGxObject go = null;
            openRaterDLg.DoModalOpen(this.Handle.ToInt32(), out go);
            if (go == null)
                return;
            go.Reset();
            IGxObject gxObject = go.Next();
            while (gxObject != null)
            {
                IGxDataset dataset = gxObject as IGxDataset;
                IRasterLayer layer = new RasterLayerClass();
                layer.CreateFromDataset(dataset.Dataset as IRasterDataset);
                this.axMapControl1.AddLayer(layer, 0);
                gxObject = go.Next();
            }

            this.axMapControl1.Refresh();
        }

        private void SelectRarsterLayerToolStripMenuItem_Click(object sender, EventArgs e)
        {
            List<ILayer> layers = new List<ILayer>();
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                layers.Add(axMapControl1.get_Layer(i));
            }
            frmSelectRasterLayer sl = new frmSelectRasterLayer(layers.ToArray());
            if (sl.ShowDialog() == DialogResult.OK)
            {
                pFirstRasterLayer = sl.rasterLayer1996;
                pSecRasterLayer = sl.rasterLayer2005;
                pForecastRasterLayer = sl.forecastLayer;
                pForecastFeatureLayer = sl.pointLayer;

                ltList = LoadLandTypeDictionary(Application.StartupPath + "\\地类表.xml");
                LogisticsFactors = LoadLogisticFactors(Application.StartupPath + "\\逻辑回归系数.xml");
                int index = 0;
                foreach (int i in ltList.Keys)
                {
                    dlbmToIndex.Add(i, index);
                    IndexTodlbm.Add(index, i);
                    index++;
                }
            }
        }

        private Dictionary<LandType, List<Land>> ErgodicRasterLand(IRasterLayer layer)
        {
            Dictionary<LandType, List<Land>> all = new Dictionary<LandType, List<Land>>();
            IRasterProps raster1996 = layer.Raster as IRasterProps;

            //影像的宽度和高度
            double width = raster1996.Width;
            double height = raster1996.Height;


            //根据影像设定分辨率
            double cellWidth = raster1996.MeanCellSize().X;
            double cellHeight = raster1996.MeanCellSize().Y;

            IPnt pntSize = new PntClass();
            pntSize.SetCoords(cellWidth, cellHeight);
            IPixelBlock pixelBlock = layer.Raster.CreatePixelBlock(pntSize);
            for (int i = 0; i < height; ++i)
            {
                for (int j = 0; j < width; ++j)
                {
                    string str = "影像分辨率为:" + height + "×" + width;
                    str += ";正在处理第" + i + "行,第" + j + "列";
                    this.statusBarXY.Text = str;
                    System.Windows.Forms.Application.DoEvents();
                    int value = Function.GetPixelValue(j, i, layer, pixelBlock);
                    if (value <= 0)
                        continue;
                    LandType lt = ltList[value];
                    if (all.ContainsKey(lt))
                        all[lt].Add(new Land(j, i, lt.dlbm));
                    else
                    {
                        List<Land> landList = new List<Land>();
                        landList.Add(new Land(j, i, lt.dlbm));
                        all.Add(lt, landList);
                    }
                }
            }
            return all;
        }

        private Dictionary<LandType, List<Land>> ErgodicRasterLand(IRasterLayer layer, Dictionary<Position, Land> ll)
        {
            Dictionary<LandType, List<Land>> all = new Dictionary<LandType, List<Land>>();
            IRasterProps raster1996 = layer.Raster as IRasterProps;

            //影像的宽度和高度
            double width = raster1996.Width;
            double height = raster1996.Height;
            //根据影像设定分辨率
            double cellWidth = raster1996.MeanCellSize().X;
            double cellHeight = raster1996.MeanCellSize().Y;

            IPnt pntSize = new PntClass();
            pntSize.SetCoords(cellWidth, cellHeight);
            IPixelBlock pixelBlock = layer.Raster.CreatePixelBlock(pntSize);
            for (int i = 0; i < height; ++i)
            {
                for (int j = 0; j < width; ++j)
                {
                    string str = "影像分辨率为:" + height + "×" + width;
                    str += ";正在处理第" + i + "行,第" + j + "列";
                    this.statusBarXY.Text = str;
                    System.Windows.Forms.Application.DoEvents();
                    int value = Function.GetPixelValue(j, i, layer, pixelBlock);
                    if (value <= 0)
                        continue;
                    LandType lt = ltList[value];
                    Land l = new Land(j, i, lt.dlbm);
                    ll.Add(new Position(j, i), l);
                    if (all.ContainsKey(lt))
                    {
                        all[lt].Add(l);
                    }
                    else
                    {
                        List<Land> landList = new List<Land>();
                        landList.Add(l);
                        all.Add(lt, landList);
                    }
                }
            }
            return all;
        }
        private Dictionary<int, LandType> LoadLandTypeDictionary(string path)
        {
            Dictionary<int, LandType> ll = new Dictionary<int, LandType>();
            XmlDocument xmlDoc = new XmlDocument();
            xmlDoc.Load(path);
            foreach (XmlNode xmlnode in xmlDoc.LastChild.ChildNodes)
            {
                int code = Convert.ToInt32(xmlnode.FirstChild.InnerText);
                string name = xmlnode.LastChild.InnerText;
                ll.Add(code, new LandType(code, name));
            }
            return ll;
        }

        private Dictionary<int, Dictionary<int, double>> LoadFactors(string path)
        {
            Dictionary<int, Dictionary<int, double>> table = new Dictionary<int, Dictionary<int, double>>();
            XmlDocument xmlDoc = new XmlDocument();
            xmlDoc.Load(path);
            foreach (XmlNode xmlnode in xmlDoc.LastChild.ChildNodes)
            {
                int code1 = Convert.ToInt32(xmlnode.Name.Split('_')[1]);
                Dictionary<int, double> tt = new Dictionary<int, double>();
                table.Add(code1, tt);
                foreach (XmlNode node in xmlnode.ChildNodes)
                {
                    tt.Add(Convert.ToInt32(node.Name.Split('_')[1]), Convert.ToDouble(node.InnerText));
                }
            }

            return table;
        }

        private Dictionary<int, Dictionary<int, LogisticFactor>> LoadLogisticFactors(string path)
        {
            Dictionary<int, Dictionary<int, LogisticFactor>> factors = new Dictionary<int, Dictionary<int, LogisticFactor>>();
            XmlDocument xmlDoc = new XmlDocument();
            xmlDoc.Load(path);

            foreach (XmlNode xmlnode in xmlDoc.LastChild.ChildNodes)
            {
                int code1 = Convert.ToInt32(xmlnode.Name.Split('_')[1]);
                Dictionary<int, LogisticFactor> tt = new Dictionary<int, LogisticFactor>();
                factors.Add(code1, tt);

                foreach (XmlNode node in xmlnode.ChildNodes)
                {
                    LogisticFactor lf = new LogisticFactor();
                    int count = 0;
                    int code2 = Convert.ToInt32(node.Name.Split('_')[1]);
                    lf.landBM = code2;
                    foreach (XmlNode n in node.ChildNodes)
                    {
                        if (count == 0)
                            lf.a = Convert.ToDouble(n.InnerText);
                        else
                            lf.b.Add(Convert.ToDouble(n.InnerText));
                        count++;
                    }
                    tt.Add(code2, lf);
                }
            }

            return factors;
        }
        private void ExportLogisticXYToolStripMenuItem_Click(object sender, EventArgs e)
        {
            frmExportFactors ef = new frmExportFactors();
            //是否已经有了自变量文件，有就从文件中读取
            if (ef.ShowDialog() == DialogResult.Cancel)
                return;
            if (ef.isHaveFiles == true)
            {
                directoryPath = ef.FileNames[0].Substring(0, ef.FileNames[0].LastIndexOf("\\"));
                MessageBox.Show("设置成功。");
            }
            else
            {
                if (!Directory.Exists(Application.StartupPath + "\\logistic"))
                {
                    Directory.CreateDirectory(Application.StartupPath + "\\logistic");
                }
                else
                {
                    Directory.Delete(Application.StartupPath + "\\logistic", true);
                    Directory.CreateDirectory(Application.StartupPath + "\\logistic");
                }

                FileStream[] fs = new FileStream[ltList.Count];
                StreamWriter[] sw = new StreamWriter[ltList.Count];
                IRaster2 rst2First = pFirstRasterLayer.Raster as IRaster2;
                IRaster2 rst2Sec = pSecRasterLayer.Raster as IRaster2;
                int index = 0;

                DateTime time1 = DateTime.Now;
                foreach (int i in ltList.Keys)
                {
                    fs[index] = new FileStream(Application.StartupPath + "\\logistic\\" + i.ToString() + ".txt", FileMode.Create);
                    sw[index] = new StreamWriter(fs[index]);
                    string tileStr = "";
                    foreach (int j in ltList.Keys)
                    {
                        tileStr += "Y(" + i.ToString() + "_" + j.ToString() + "),";
                    }
                    tileStr += "Population,Elevation,Slope,Dist(商业中心),Dist(金融中心),Dist(工业中心),Dist(文教中心),Dist(Railways),Dist(Roads)";
                    sw[index].WriteLine(tileStr);
                    index++;
                }

                #region 自变量图层
                IFeatureLayer pPop = ArcEngineFunction.GetLayerByName(axMapControl1, "Population_value") as IFeatureLayer;
                IFeatureLayer pDem = ArcEngineFunction.GetLayerByName(axMapControl1, "Dem_value") as IFeatureLayer;
                IFeatureLayer pSlope = ArcEngineFunction.GetLayerByName(axMapControl1, "Slope_value") as IFeatureLayer;
                IFeatureLayer pSY = ArcEngineFunction.GetLayerByName(axMapControl1, "shopping_dist") as IFeatureLayer;
                IFeatureLayer pJR = ArcEngineFunction.GetLayerByName(axMapControl1, "金融中心_dist") as IFeatureLayer;
                IFeatureLayer pGY = ArcEngineFunction.GetLayerByName(axMapControl1, "工业中心_dist") as IFeatureLayer;
                IFeatureLayer pWJ = ArcEngineFunction.GetLayerByName(axMapControl1, "school_dist") as IFeatureLayer;
                IFeatureLayer pRailway = ArcEngineFunction.GetLayerByName(axMapControl1, "Railway_dist") as IFeatureLayer;
                IFeatureLayer pRoad = ArcEngineFunction.GetLayerByName(axMapControl1, "0-6级道路中心线_dist") as IFeatureLayer;
                #endregion
                IRasterProps raster1996 = pForecastRasterLayer.Raster as IRasterProps;
                //影像的宽度和高度
                double width = raster1996.Width;
                double height = raster1996.Height;
                int cellcounts = pForecastFeatureLayer.FeatureClass.FeatureCount(null);
                cellValue = new int[cellcounts];
                for (int i = 0; i < cellcounts; i++)
                {

                    IFeature f = pForecastFeatureLayer.FeatureClass.GetFeature(i);
                    IPoint p = f.Shape as IPoint;
                    int row = rst2First.ToPixelRow(p.Y);
                    int col = rst2First.ToPixelColumn(p.X);
                    int value1 = Convert.ToInt32(rst2First.GetPixelValue(0, col, row));
                    int ind = dlbmToIndex[value1];
                    int value2 = Convert.ToInt32(rst2Sec.GetPixelValue(0, col, row));

                    string str1 = "影像共有:" + cellcounts.ToString() + "个元胞";
                    str1 += ";正在处理第" + i.ToString() + "个。";
                    this.statusBarXY.Text = str1;
                    System.Windows.Forms.Application.DoEvents();
                    //Land l = new Land(col, row, value1);

                    string str = "";
                    foreach (int code in ltList.Keys)
                    {
                        if (value2 == code)
                        {
                            str += "1,";
                        }
                        else
                        {
                            str += "0,";
                        }
                    }
                    string factors = "";

                    #region 自变量
                    //从字段中读取
                    double filedValue = Convert.ToDouble(pPop.FeatureClass.GetFeature(i).get_Value(4));
                    factors += filedValue.ToString() + ",";
                    //l.XFactors.Add(filedValue);

                    filedValue = Convert.ToDouble(pDem.FeatureClass.GetFeature(i).get_Value(4));
                    factors += filedValue.ToString() + ",";
                    //l.XFactors.Add(filedValue);

                    filedValue = Convert.ToDouble(pSlope.FeatureClass.GetFeature(i).get_Value(4));
                    factors += filedValue.ToString() + ",";
                    //l.XFactors.Add(filedValue);

                    filedValue = Convert.ToDouble(pSY.FeatureClass.GetFeature(i).get_Value(5));
                    factors += filedValue.ToString() + ",";
                    //l.XFactors.Add(filedValue);

                    filedValue = Convert.ToDouble(pJR.FeatureClass.GetFeature(i).get_Value(5));
                    factors += filedValue.ToString() + ",";
                    //l.XFactors.Add(filedValue);

                    filedValue = Convert.ToDouble(pGY.FeatureClass.GetFeature(i).get_Value(5));
                    factors += filedValue.ToString() + ",";
                    //l.XFactors.Add(filedValue);

                    filedValue = Convert.ToDouble(pWJ.FeatureClass.GetFeature(i).get_Value(5));
                    factors += filedValue.ToString() + ",";
                    //l.XFactors.Add(filedValue);

                    filedValue = Convert.ToDouble(pRailway.FeatureClass.GetFeature(i).get_Value(5));
                    factors += filedValue.ToString() + ",";
                    //l.XFactors.Add(filedValue);

                    filedValue = Convert.ToDouble(pRoad.FeatureClass.GetFeature(i).get_Value(5));
                    factors += filedValue.ToString();
                    //l.XFactors.Add(filedValue);

                    #endregion

                    //posLand.Add(new Position(row,col), l);
                    str += factors;
                    sw[ind].WriteLine(str);


                }
                index = 0;
                foreach (int i in ltList.Keys)
                {
                    sw[index].Flush();
                    sw[index].Close();
                    fs[index].Close();
                    index++;
                }

                DateTime time2 = DateTime.Now;
                MessageBox.Show((time2 - time1).ToString());
                MessageBox.Show("导出完成。");
            }

        }

        private void ImportLogisticFactorToolStripMenuItem_Click(object sender, EventArgs e)
        {

            #region 适应度
            //foreach (LandType lt in forecastAllLand.Keys)
            //{
            //    foreach (Land l in forecastAllLand[lt])
            //    {

            //    }
            //}
            #endregion
        }

        private void CalcTotalDemandToolStripMenuItem_Click(object sender, EventArgs e)
        {
            Dictionary<int, Dictionary<int, double>> NeiborFactors = new Dictionary<int, Dictionary<int, double>>();//邻居偏好系数表
            Dictionary<int, Dictionary<int, double>> InheritFactors = new Dictionary<int, Dictionary<int, double>>();//遗传系数表
            NeiborFactors = LoadFactors(Application.StartupPath + "\\偏好系数.xml");
            InheritFactors = LoadFactors(Application.StartupPath + "\\遗传系数.xml");
            IFeatureLayer pSTX = ArcEngineFunction.GetLayerByName(axMapControl1, "生态控制线") as IFeatureLayer;//生态线图层
            IFeature stx = pSTX.FeatureClass.GetFeature(0);
            IRelationalOperator ro = stx.Shape as IRelationalOperator;
            int cellcounts = pForecastFeatureLayer.FeatureClass.FeatureCount(null);
            bool isBreak = true;
            int num = 1;
            Dictionary<string, int> AllCellValues = new Dictionary<string, int>();
            Dictionary<string, int> AllCellValuesCopy = new Dictionary<string, int>();//备份，使在计算时不改变元胞的值
            Dictionary<string, string> AllFactors = new Dictionary<string, string>();
            int CycleCount = 0;
            #region 读取像元值和自变量因数
            string filepath = directoryPath + "\\Y-X.txt";
            FileStream fs = new FileStream(filepath, FileMode.Open);
            StreamReader sr = new StreamReader(fs);
            sr.ReadLine();
            IRasterProps tmpRaster = pForecastRasterLayer.Raster as IRasterProps;
            IRaster2 rst2 = pForecastRasterLayer.Raster as IRaster2;
            for (int i = 0; i < cellcounts; i++)
            {
                //以12年为模拟基础，因为12年比09年少一个象元，所以在这个要素这要跳过，如果是以09年为模拟基础就不用这段
                if (i == 2095414)
                {
                    sr.ReadLine();
                    continue;
                }

                //string str1 = "正在读取像元，共有:" + cellcounts.ToString() + "个元胞";
                //str1 += ";正在处理第" + (i + 1) + "个。";
                //this.statusBarXY.Text = str1;
                //System.Windows.Forms.Application.DoEvents();
                IFeature f = pForecastFeatureLayer.FeatureClass.GetFeature(i);
                IPoint p = f.Shape as IPoint;
                int row = rst2.ToPixelRow(p.Y);
                int col = rst2.ToPixelColumn(p.X);
                int value1 = Convert.ToInt32(rst2.GetPixelValue(0, col, row));
                string rowcol = row + "_" + col;
                AllCellValues.Add(rowcol, value1);
                AllCellValuesCopy.Add(rowcol, value1);
                string str = sr.ReadLine();
                AllFactors.Add(rowcol, str);
            }
            fs.Close();
            sr.Close();
            #endregion
            while (isBreak)
            {
                CycleCount++;
                //int BreakCount = 0;//中断条件，用地类型221的数量
                //X变量文件
                Dictionary<string, double[]> AllFitness = new Dictionary<string, double[]>();
                int index = 0;
                IRaster2 rst2First = pForecastRasterLayer.Raster as IRaster2;
                IRaster raster = pForecastRasterLayer.Raster;
                IRasterProps raster1996 = raster as IRasterProps;

                //   指定像素块大小
                IPnt pntBlock = new PntClass();
                pntBlock.X = raster1996.Width;
                pntBlock.Y = raster1996.Height;
                //创建一个光标以给定像素块大小
                IRasterCursor rasterCursor = rst2First.CreateCursorEx(pntBlock);
                //创建像素块
                IPixelBlock3 pixelBlock3 = null;
                int blockWidth = 0;
                int blockHeight = 0;
                pixelBlock3 = rasterCursor.PixelBlock as IPixelBlock3;
                blockWidth = pixelBlock3.Width;
                blockHeight = pixelBlock3.Height;
                int width = tmpRaster.Width;
                int height = tmpRaster.Height;
                List<string> ChangeRowCol = new List<string>();//记录已经发生改变了的行列号
                //指定平面的像素的数组
                System.Array pixels = (System.Array)pixelBlock3.get_PixelData(0);
                for (int i = 0; i < cellcounts; i++)
                {
                    //string str1 = "影像共有:" + cellcounts + "个元胞" + ";正在处理第" + (i + 1) + "个。";
                    //this.statusBarXY.Text = str1;
                    //System.Windows.Forms.Application.DoEvents();
                    if (i == 2095414)
                        continue;
                    IFeature f = pForecastFeatureLayer.FeatureClass.GetFeature(i);
                    IPoint p = f.Shape as IPoint;
                    int row = rst2First.ToPixelRow(p.Y);
                    int col = rst2First.ToPixelColumn(p.X);
                    string rowcol = row + "_" + col;
                    int value1 = AllCellValues[rowcol];
                    double[] fits = GetFitness(rowcol, AllCellValues, AllFactors, AllFitness);

                    double[] totalDemand = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
                    double[] NeiDemand = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
                    double[] InheritDemand = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };

                    int index1 = 0;
                    foreach (int j in ltList.Keys)
                    {
                        InheritDemand[index1] = fits[index1] * InheritFactors[value1][j];
                        index1++;
                    }

                    #region 邻居欲望
                    int neiCount = 0;
                    string tmprowcol = (row - 1) + "_" + (col - 1);
                    if (AllCellValues.ContainsKey(tmprowcol))
                    {
                        neiCount++;
                        index = 0;
                        int value = AllCellValues[tmprowcol];
                        double[] fitness = GetFitness(tmprowcol, AllCellValues, AllFactors, AllFitness);
                        foreach (int j in ltList.Keys)
                        {
                            NeiDemand[index] += fitness[index] * NeiborFactors[value][j];
                            index++;
                        }
                    }
                    //上
                    //pos = new Position(l.x , l.y - 1);
                    tmprowcol = (row - 1) + "_" + (col);
                    if (AllCellValues.ContainsKey(tmprowcol))
                    {
                        neiCount++;
                        index = 0;
                        int value = AllCellValues[tmprowcol];
                        double[] fitness = GetFitness(tmprowcol, AllCellValues, AllFactors, AllFitness);
                        foreach (int j in ltList.Keys)
                        {
                            NeiDemand[index] += fitness[index] * NeiborFactors[value][j];
                            index++;
                        }
                    }

                    //右上
                    //pos = new Position(l.x+1, l.y - 1);
                    tmprowcol = (row - 1) + "_" + (col + 1);
                    if (AllCellValues.ContainsKey(tmprowcol))
                    {
                        neiCount++;
                        index = 0;
                        int value = AllCellValues[tmprowcol];
                        double[] fitness = GetFitness(tmprowcol, AllCellValues, AllFactors, AllFitness);
                        foreach (int j in ltList.Keys)
                        {
                            NeiDemand[index] += fitness[index] * NeiborFactors[value][j];
                            index++;
                        }
                    }

                    //右
                    //pos = new Position(l.x+1, l.y);
                    tmprowcol = (row) + "_" + (col + 1);
                    if (AllCellValues.ContainsKey(tmprowcol))
                    {
                        neiCount++;
                        index = 0;
                        int value = AllCellValues[tmprowcol];
                        double[] fitness = GetFitness(tmprowcol, AllCellValues, AllFactors, AllFitness);
                        foreach (int j in ltList.Keys)
                        {
                            NeiDemand[index] += fitness[index] * NeiborFactors[value][j];
                            index++;
                        }
                    }



                    //右下
                    //pos = new Position(l.x+1, l.y +1);
                    tmprowcol = (row + 1) + "_" + (col + 1);
                    if (AllCellValues.ContainsKey(tmprowcol))
                    {
                        neiCount++;
                        index = 0;
                        int value = AllCellValues[tmprowcol];
                        double[] fitness = GetFitness(tmprowcol, AllCellValues, AllFactors, AllFitness);
                        foreach (int j in ltList.Keys)
                        {
                            NeiDemand[index] += fitness[index] * NeiborFactors[value][j];
                            index++;
                        }
                    }


                    //下
                    //pos = new Position(l.x, l.y + 1);
                    tmprowcol = (row + 1) + "_" + (col);
                    if (AllCellValues.ContainsKey(tmprowcol))
                    {
                        neiCount++;
                        index = 0;
                        int value = AllCellValues[tmprowcol];
                        double[] fitness = GetFitness(tmprowcol, AllCellValues, AllFactors, AllFitness);
                        foreach (int j in ltList.Keys)
                        {
                            NeiDemand[index] += fitness[index] * NeiborFactors[value][j];
                            index++;
                        }
                    }



                    //左下
                    //pos = new Position(l.x-1, l.y+1);
                    tmprowcol = (row + 1) + "_" + (col - 1);
                    if (AllCellValues.ContainsKey(tmprowcol))
                    {
                        neiCount++;
                        index = 0;
                        int value = AllCellValues[tmprowcol];
                        double[] fitness = GetFitness(tmprowcol, AllCellValues, AllFactors, AllFitness);
                        foreach (int j in ltList.Keys)
                        {
                            NeiDemand[index] += fitness[index] * NeiborFactors[value][j];
                            index++;
                        }
                    }



                    //左 
                    //pos = new Position(l.x - 1, l.y);
                    tmprowcol = (row) + "_" + (col - 1);
                    if (AllCellValues.ContainsKey(tmprowcol))
                    {
                        neiCount++;
                        index = 0;
                        int value = AllCellValues[tmprowcol];
                        double[] fitness = GetFitness(tmprowcol, AllCellValues, AllFactors, AllFitness);
                        foreach (int j in ltList.Keys)
                        {
                            NeiDemand[index] += fitness[index] * NeiborFactors[value][j];
                            index++;
                        }
                    }

                    #endregion


                    for (int j = 0; j < totalDemand.Length; j++)
                    {
                        totalDemand[j] = (Math.Pow(1.03, CycleCount) * NeiDemand[j] / neiCount) + (Math.Pow(0.99, CycleCount) * InheritDemand[j]);
                    }

                    int tmpIndex = FindMaxIndex(totalDemand);
                    int changevalue = IndexTodlbm[tmpIndex - 1];
                    //AllCellValues[rowcol]= changevalue;

                    #region 判断元胞是否转换
                    bool isChange = false;
                    if (changevalue != value1)
                    {
                        if (fjsyd.Contains(value1))
                        {
                            if (fjsyd.Contains(changevalue))
                            {
                                //if (value1 != changevalue)
                                isChange = true;
                            }
                            else
                            {
                                if (!IsContains(p, pSTX as IFeatureLayer))
                                {
                                    //if (value1 != changevalue)
                                    isChange = true;
                                }
                            }
                        }
                        else
                        {
                            //if (value1 != changevalue)
                            isChange = true;
                        }
                    }
                    #endregion
                    if (isChange)
                    {
                        AllCellValuesCopy[rowcol] = changevalue;
                        pixels.SetValue(Convert.ToByte(changevalue), col, row);
                        ChangeRowCol.Add(rowcol);
                    }
                    else
                    {
                        pixels.SetValue(Convert.ToByte(value1), col, row);
                    }
                }

                IRasterEdit rasterEdit = rst2First as IRasterEdit;
                if (rasterEdit.CanEdit())
                {

                    pixelBlock3.set_PixelData(0, pixels);
                    IPnt pnt = null;
                    pnt = rasterCursor.TopLeft;
                    rasterEdit.Write(pnt, (IPixelBlock)pixelBlock3);
                    rasterEdit.Refresh();
                    IRaster pRaster = pForecastRasterLayer.Raster;
                    IRaster2 pRaster2 = pRaster as IRaster2;
                    ISaveAs pSaveAs = pRaster2 as ISaveAs;
                    pSaveAs.SaveAs("E:\\s私活\\AETools\\Data\\栅格备份\\Raster2009_" + num + ".tif", null, "TIFF");
                    num++;
                }

                foreach (string item in ChangeRowCol)
                {
                    AllCellValues[item] = AllCellValuesCopy[item];
                }
                ChangeRowCol.Clear();
                System.Runtime.InteropServices.Marshal.ReleaseComObject(rasterEdit);

                //if (Math.Abs(BreakCount - 313109) / 313109 <= 0.005)
                //{
                //    isBreak = true;
                //    IRaster pRaster = pForecastRasterLayer.Raster;
                //    IRaster2 pRaster2 = pRaster as IRaster2;
                //    ISaveAs pSaveAs = pRaster2 as ISaveAs;
                //    pSaveAs.SaveAs("E:\\s私活\\AETools\\Data\\栅格备份\\result.tif", null, "TIFF");
                //}
            }
            //MessageBox.Show("模拟完成！");
        }

        /// <summary>
        /// 计算地块适应度
        /// </summary>
        /// <param name="rowcol">行列号</param>
        /// <param name="cvs">cellvalues</param>
        /// <param name="fs">factors</param>
        /// <param name="fits">fitness</param>
        /// <returns></returns>
        private double[] GetFitness(string rowcol, Dictionary<string, int> cvs, Dictionary<string, string> fs, Dictionary<string, double[]> fits)
        {
            string str = fs[rowcol];
            int value1 = cvs[rowcol];
            string[] factors = str.Split(',');
            //if (fits.ContainsKey(rowcol))
            //    return fits[rowcol];
            //else
            //{
            double[] fitness = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };//地块适应度
            int ind = 0;
            foreach (int j in ltList.Keys)
            {
                LogisticFactor lf = LogisticsFactors[value1][j];
                fitness[ind] = lf.a;
                for (int k = 0; k < lf.b.Count; k++)
                {
                    fitness[ind] += lf.b[k] * Convert.ToDouble(factors[k + 16]);
                }
                if (fitness[ind] != 0)
                    fitness[ind] = 1 / (1 + Math.Exp(-fitness[ind]));
                ind++;
            }
            //fits.Add(rowcol, fitness);
            return fitness;
            //}

        }


        private bool IsInRasterExt(int col, int row, int allrow, int allcol)
        {
            if (row >= 0 && row < allrow & col >= 0 && col < allcol)
                return true;
            else
                return false;
        }

        private double GetShortestDistance(IGeometry g, IFeatureLayer layer)
        {
            double dist = 0;

            IFeatureCursor pCursor = layer.Search(null, false);
            IFeature pFeature;
            pFeature = pCursor.NextFeature();
            IProximityOperator pProximity = (IProximityOperator)g;
            dist = pProximity.ReturnDistance(pFeature.Shape);
            int count = 0;
            while (pFeature != null)
            {
                count++;
                double dd = pProximity.ReturnDistance(pFeature.Shape);
                if (dist >= dd)
                    dist = dd;
                pFeature = pCursor.NextFeature();
            }
            return dist;
        }

        private double GetShortestDistance(IGeometry g, IFeatureLayer layer, bool one)
        {
            double dist = 0;
            IProximityOperator pProximity = (IProximityOperator)g;

            dist = pProximity.ReturnDistance(layer.FeatureClass.GetFeature(0).Shape);
            return dist;
        }

        private bool IsContains(IGeometry g, IFeatureLayer layer)
        {
            IFeatureCursor pCursor = layer.Search(null, false);
            IFeature pFeature;
            pFeature = pCursor.NextFeature();
            IRelationalOperator pProximity = (IRelationalOperator)g;
            while (pFeature != null)
            {
                if (pProximity.Contains(pFeature.Shape))
                    return true;
                pFeature = pCursor.NextFeature();
            }
            return false;
        }
        private void SelectDataAndExportResultToolStripMenuItem_Click(object sender, EventArgs e)
        {
        }

        private int FindMaxIndex(double[] a)
        {
            int index = 0;
            double value = a[0];
            for (int i = 1; i < a.Length; i++)
            {
                if (value < a[i])
                {
                    value = a[i];
                    index = i;
                }
            }
            return (index + 1);
        }

        private void 马尔科夫概率计算ToolStripMenuItem_Click(object sender, EventArgs e)
        {

            List<ILayer> layers = new List<ILayer>();
            IRasterLayer layer1 = null;//第一次影像
            IRasterLayer layer2 = null;//第二次影像

            Dictionary<int, Dictionary<int, int>> Markov = new Dictionary<int, Dictionary<int, int>>();
            Dictionary<int, int> landCount = new Dictionary<int, int>();//统计地块总数
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                layers.Add(axMapControl1.get_Layer(i));
            }
            frmSelectRasterLayer sl = new frmSelectRasterLayer(layers.ToArray());
            if (sl.ShowDialog() == DialogResult.OK)
            {
                layer1 = sl.rasterLayer1996;
                layer2 = sl.rasterLayer2005;

                IRasterProps raster1996 = layer1.Raster as IRasterProps;
                IRasterProps raster2005 = layer2.Raster as IRasterProps;

                double width = raster1996.Width;
                double height = raster1996.Height;


                //根据影像设定分辨率
                double cellWidth = raster1996.MeanCellSize().X;
                double cellHeight = raster1996.MeanCellSize().Y;

                IPnt pntSize = new PntClass();
                pntSize.SetCoords(cellWidth, cellHeight);
                IPixelBlock pixelBlock = layer1.Raster.CreatePixelBlock(pntSize);
                for (int i = 0; i < height; ++i)
                {
                    //DateTime dt1 = DateTime.Now;
                    for (int j = 0; j < width; ++j)
                    {
                        string str = "影像分辨率为:" + height + "×" + width;
                        str += ";正在处理第" + i + "行,第" + j + "列";
                        this.statusBarXY.Text = str;
                        System.Windows.Forms.Application.DoEvents();
                        int value1 = Function.GetPixelValue(j, i, layer1, pixelBlock);
                        if (value1 <= 0)
                            continue;
                        int value2 = Function.GetPixelValue(j, i, layer2, pixelBlock);
                        if (value2 <= 0)
                            continue;
                        if (landCount.ContainsKey(value1))
                        {
                            landCount[value1]++;
                        }
                        else
                        {
                            landCount.Add(value1, 1);
                        }
                        if (Markov.ContainsKey(value1))
                        {
                            Dictionary<int, int> tmpDic = new Dictionary<int, int>();
                            tmpDic = Markov[value1];
                            if (tmpDic.ContainsKey(value2))
                            {
                                tmpDic[value2]++;
                            }
                            else
                            {
                                tmpDic.Add(value2, 1);
                            }
                        }
                        else
                        {
                            Dictionary<int, int> tmpDic = new Dictionary<int, int>();
                            tmpDic.Add(value2, 1);
                            Markov.Add(value1, tmpDic);
                        }

                    }
                    //DateTime dt2 = DateTime.Now;
                    //MessageBox.Show((dt2 - dt1).ToString());
                }

                XmlDocument xmlDoc = new XmlDocument();
                XmlElement root = xmlDoc.CreateElement("Markov"); // 创建根节点album
                xmlDoc.AppendChild(root);
                foreach (int dlbm1 in Markov.Keys)
                {
                    XmlElement xmlA = xmlDoc.CreateElement("_" + dlbm1.ToString()); // 创建根节点album
                    root.AppendChild(xmlA);
                    foreach (int dlbm2 in Markov[dlbm1].Keys)
                    {
                        XmlElement xmlB = xmlDoc.CreateElement("_" + dlbm2.ToString()); // 创建根节点album
                        xmlA.AppendChild(xmlB);
                        double aa = (double)Markov[dlbm1][dlbm2] / (double)landCount[dlbm1];
                        aa = Math.Round(aa, 4);
                        xmlB.InnerText = aa.ToString();
                    }
                }

                xmlDoc.Save(Application.StartupPath + "\\markov.xml");
                MessageBox.Show("计算完成。");
            }
        }

        private void 栅格中心元胞比例ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            List<ILayer> layers = new List<ILayer>();
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                layers.Add(axMapControl1.get_Layer(i));
            }
            frmSelectRasterLayer sl = new frmSelectRasterLayer(layers.ToArray());

            if (sl.ShowDialog() == DialogResult.OK)
            {
                IRasterLayer pRL1 = null;
                IRasterLayer pRL2 = null;
                pRL1 = sl.rasterLayer1996;
                pRL2 = sl.rasterLayer2005;
                IRasterProps pRP1 = pRL1.Raster as IRasterProps;
                IRasterProps pRP2 = pRL2.Raster as IRasterProps;
                //比较两年的栅格数据的分辨率是否一致
                //不一致就返回
                if (pRP1.MeanCellSize().X != pRP2.MeanCellSize().X ||
                    pRP1.MeanCellSize().Y != pRP2.MeanCellSize().Y)
                {
                    string str = "所选两年的栅格数据的分辨率不一致,请重新选择！\r\n";
                    MessageBox.Show(str);
                    return;
                }
                ArcEngineFunction.ResetAllByCalc11X11(pRL1);
                ArcEngineFunction.ResetAllByCalc11X11(pRL2);
                MessageBox.Show("计算结束。");
            }
        }

        private void 提取扩展部分的栅格ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            List<ILayer> layers = new List<ILayer>();
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                layers.Add(axMapControl1.get_Layer(i));
            }
            frmSelectRasterLayer sl = new frmSelectRasterLayer(layers.ToArray());

            if (sl.ShowDialog() == DialogResult.OK)
            {
                IRasterLayer pRL1 = null;
                IRasterLayer pRL2 = null;
                pRL1 = sl.rasterLayer1996;
                pRL2 = sl.rasterLayer2005;
                IRasterProps pRP1 = pRL1.Raster as IRasterProps;
                IRasterProps pRP2 = pRL2.Raster as IRasterProps;
                //比较两年的栅格数据的分辨率是否一致
                //不一致就返回
                if (pRP1.MeanCellSize().X != pRP2.MeanCellSize().X ||
                    pRP1.MeanCellSize().Y != pRP2.MeanCellSize().Y)
                {
                    string str = "所选两年的栅格数据的分辨率不一致,请重新选择！\r\n";
                    MessageBox.Show(str);
                    return;
                }
                pRP1.PixelType = rstPixelType.PT_DOUBLE;
                //影像的宽度和高度
                double width = pRP1.Width;
                double height = pRP1.Height;


                //根据影像设定分辨率
                double cellWidth = pRP1.MeanCellSize().X;
                double cellHeight = pRP1.MeanCellSize().Y;
                IPnt pntSize = new PntClass();
                pntSize.SetCoords(cellWidth, cellHeight);
                IPixelBlock pixelBlock = pRL2.Raster.CreatePixelBlock(pntSize);
                double cellSize = cellWidth;
                IRaster2 rst2First = pRL2.Raster as IRaster2;
                IRaster raster = pRL2.Raster;
                IRasterProps raster1996 = raster as IRasterProps;

                //   指定像素块大小
                IPnt pntBlock = new PntClass();
                pntBlock.X = raster1996.Width;
                pntBlock.Y = raster1996.Height;
                //创建一个光标以给定像素块大小
                IRasterCursor rasterCursor = rst2First.CreateCursorEx(pntBlock);
                //创建像素块
                IPixelBlock3 pixelBlock3 = null;
                int blockWidth = 0;
                int blockHeight = 0;
                pixelBlock3 = rasterCursor.PixelBlock as IPixelBlock3;
                blockWidth = pixelBlock3.Width;
                blockHeight = pixelBlock3.Height;
                double averange = 0;
                int count = 0;
                //指定平面的像素的数组
                System.Array pixels = (System.Array)pixelBlock3.get_PixelData(0);
                for (int i = 0; i < height; ++i)
                {
                    for (int j = 0; j < width; ++j)
                    {
                        //double value2 = ArcEngineFunction.GetPxValue(j, i, pRL2, pixelBlock);
                        IPnt pnt = new PntClass();
                        pnt.SetCoords(j, i);
                        pRL2.Raster.Read(pnt, pixelBlock);
                        object obj = pixelBlock.GetVal(0, 0, 0);
                        if (obj == null)
                            continue;
                        count++;
                        averange += Convert.ToDouble(obj);
                    }
                }
                averange /= count;
                //double nodatavalue=-1.797693*Math.Pow(10,38);
                //averange = 0.123925;
                for (int i = 0; i < height; ++i)
                {
                    for (int j = 0; j < width; ++j)
                    {
                        //double value1 = ArcEngineFunction.GetPxValue(j, i, pRL1, pixelBlock);
                        //double value2 = ArcEngineFunction.GetPxValue(j, i, pRL2, pixelBlock);
                        IPnt pnt = new PntClass();
                        pnt.SetCoords(j, i);
                        pRL1.Raster.Read(pnt, pixelBlock);
                        object obj1 = pixelBlock.GetVal(0, 0, 0);
                        pRL2.Raster.Read(pnt, pixelBlock);
                        object obj2 = pixelBlock.GetVal(0, 0, 0);
                        if (obj1 == null && obj2 != null)
                        {
                            double result = Convert.ToDouble(obj2) / averange;
                            result = Math.Round(result, 2);
                            if (result >= 1)
                                pixels.SetValue(3.0, j, i);
                            else
                                pixels.SetValue(2.0, j, i);
                        }
                        else
                        {
                            //-1.797693e+308
                            pixels.SetValue(0.0, j, i);
                        }

                    }
                }

                IRasterEdit rasterEdit = rst2First as IRasterEdit;
                if (rasterEdit.CanEdit())
                {

                    pixelBlock3.set_PixelData(0, pixels);
                    IPnt pnt = null;
                    pnt = rasterCursor.TopLeft;
                    rasterEdit.Write(pnt, (IPixelBlock)pixelBlock3);
                    rasterEdit.Refresh();
                }

                System.Runtime.InteropServices.Marshal.ReleaseComObject(rasterEdit);


            }
            MessageBox.Show("处理完成。");
        }

        private void 差值结果图ToolStripMenuItem_Click(object sender, EventArgs e)
        {

            List<ILayer> layers = new List<ILayer>();
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                layers.Add(axMapControl1.get_Layer(i));
            }
            frmSelectRasterLayer sl = new frmSelectRasterLayer(layers.ToArray());

            if (sl.ShowDialog() == DialogResult.OK)
            {
                IRasterLayer pRL1 = null;
                IRasterLayer pRL2 = null;
                pRL1 = sl.rasterLayer1996;
                pRL2 = sl.rasterLayer2005;
                IRasterProps pRP1 = pRL1.Raster as IRasterProps;
                IRasterProps pRP2 = pRL2.Raster as IRasterProps;
                //比较两年的栅格数据的分辨率是否一致
                //不一致就返回
                if (pRP1.MeanCellSize().X != pRP2.MeanCellSize().X ||
                    pRP1.MeanCellSize().Y != pRP2.MeanCellSize().Y)
                {
                    string str = "所选两年的栅格数据的分辨率不一致,请重新选择！\r\n";
                    MessageBox.Show(str);
                    return;
                }
                pRP1.PixelType = rstPixelType.PT_DOUBLE;
                //影像的宽度和高度
                double width = pRP1.Width;
                double height = pRP1.Height;


                //根据影像设定分辨率
                double cellWidth = pRP1.MeanCellSize().X;
                double cellHeight = pRP1.MeanCellSize().Y;
                IPnt pntSize = new PntClass();
                pntSize.SetCoords(cellWidth, cellHeight);
                IPixelBlock pixelBlock = pRL2.Raster.CreatePixelBlock(pntSize);
                double cellSize = cellWidth;
                IRaster2 rst2First = pRL1.Raster as IRaster2;
                IRaster raster = pRL1.Raster;
                IRasterProps raster1996 = raster as IRasterProps;

                //   指定像素块大小
                IPnt pntBlock = new PntClass();
                pntBlock.X = raster1996.Width;
                pntBlock.Y = raster1996.Height;
                //创建一个光标以给定像素块大小
                IRasterCursor rasterCursor = rst2First.CreateCursorEx(pntBlock);
                //创建像素块
                IPixelBlock3 pixelBlock3 = null;
                int blockWidth = 0;
                int blockHeight = 0;
                pixelBlock3 = rasterCursor.PixelBlock as IPixelBlock3;
                blockWidth = pixelBlock3.Width;
                blockHeight = pixelBlock3.Height;
                //指定平面的像素的数组
                System.Array pixels = (System.Array)pixelBlock3.get_PixelData(0);
                for (int i = 0; i < height; ++i)
                {
                    for (int j = 0; j < width; ++j)
                    {
                        double value1 = ArcEngineFunction.GetPxValue(j, i, pRL1, pixelBlock);
                        if (value1 == 0)
                            continue;
                        double value2 = ArcEngineFunction.GetPxValue(j, i, pRL2, pixelBlock);
                        //if (value1 == 0 && value2 != 0)
                        //{
                        //    double result = value2 / averange;
                        //    result = Math.Round(result, 2);
                        //    if (result >= 1)
                        //        pixels.SetValue(3.0, j, i);
                        //    else
                        //        pixels.SetValue(2.0, j, i);
                        //}
                        //else
                        //{
                        //    pixels.SetValue(0, j, i);
                        //}
                        double result = value2 - value1;
                        if (result <= -0.25)
                            pixels.SetValue(1.0, j, i);
                        else if (result > -0.25 && result < 0)
                            pixels.SetValue(2.0, j, i);
                        else if (result == 0)
                            pixels.SetValue(3.0, j, i);
                        else if (result < 0.25 && result > 0)
                            pixels.SetValue(4.0, j, i);
                        else if (result >= 0.25)
                            pixels.SetValue(5.0, j, i);
                        //pixels.SetValue(value2-value1, j, i);
                    }
                }

                IRasterEdit rasterEdit = rst2First as IRasterEdit;
                if (rasterEdit.CanEdit())
                {

                    pixelBlock3.set_PixelData(0, pixels);
                    IPnt pnt = null;
                    pnt = rasterCursor.TopLeft;
                    rasterEdit.Write(pnt, (IPixelBlock)pixelBlock3);
                    rasterEdit.Refresh();
                }

                System.Runtime.InteropServices.Marshal.ReleaseComObject(rasterEdit);


            }
            MessageBox.Show("处理完成。");
        }

        private void 将多个YX文件合并成一个ToolStripMenuItem_Click(object sender, EventArgs e)
        {

            frmExportFactors ef = new frmExportFactors();
            if (ef.ShowDialog() == DialogResult.Cancel)
                return;
            if (ef.isHaveFiles == true)
            {
                directoryPath = ef.FileNames[0].Substring(0, ef.FileNames[0].LastIndexOf("\\"));
                MessageBox.Show("设置成功。");
            }
            int index = 0;
            FileStream[] fs = new FileStream[16];
            StreamReader[] sr = new StreamReader[16];
            FileStream outfs = new FileStream(directoryPath + "\\Y-X.txt", FileMode.Create);
            StreamWriter outsw = new StreamWriter(outfs);
            foreach (int j in ltList.Keys)
            {
                string filepath = directoryPath + "\\" + j.ToString() + ".txt";
                fs[index] = new FileStream(filepath, FileMode.Open);
                sr[index] = new StreamReader(fs[index]);
                sr[index].ReadLine();
                index++;
            }

            IRaster2 rst2First = pForecastRasterLayer.Raster as IRaster2;
            //IProximityOperator po = stx.Shape as IRelationalOperator;


            IRaster raster = pForecastRasterLayer.Raster;
            IRasterProps raster1996 = raster as IRasterProps;

            //   指定像素块大小
            IPnt pntBlock = new PntClass();
            pntBlock.X = raster1996.Width;
            pntBlock.Y = raster1996.Height;
            //创建一个光标以给定像素块大小

            int cellcounts = pForecastFeatureLayer.FeatureClass.FeatureCount(null);
            for (int i = 0; i < cellcounts; i++)
            {
                string str1 = "影像共有:" + cellcounts.ToString() + "个元胞";
                str1 += ";正在处理第" + (i + 1) + "个。";
                this.statusBarXY.Text = str1;
                System.Windows.Forms.Application.DoEvents();
                IFeature f = pForecastFeatureLayer.FeatureClass.GetFeature(i);
                IPoint p = f.Shape as IPoint;
                int row = rst2First.ToPixelRow(p.Y);
                int col = rst2First.ToPixelColumn(p.X);
                int value1 = 0;
                value1 = Convert.ToInt32(rst2First.GetPixelValue(0, col, row));

                //读取地块参数
                int ind = 0;
                string str = "";
                ind = dlbmToIndex[value1];
                str = sr[ind].ReadLine();
                outsw.WriteLine(str);
            }
            index = 0;
            foreach (int j in ltList.Keys)
            {

                fs[index].Close();
                sr[index].Close();
                index++;
            }
            outsw.Close();
            outfs.Close();

        }


        private void toolStripMenuItem1_Click(object sender, EventArgs e)
        {

            List<ILayer> layers = new List<ILayer>();
            IRasterLayer layer1 = null;//第一次影像
            IRasterLayer layer2 = null;//第二次影像
            Dictionary<int, int> RightCount = new Dictionary<int, int>();//统计地块总数
            Dictionary<int, int> landCount = new Dictionary<int, int>();//统计地块总数
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                layers.Add(axMapControl1.get_Layer(i));
            }
            frmSelectRasterLayer sl = new frmSelectRasterLayer("预测图层", "实际图层", layers.ToArray());
            if (sl.ShowDialog() == DialogResult.OK)
            {
                layer2 = sl.rasterLayer1996;
                layer1 = sl.rasterLayer2005;

                IRasterProps raster1996 = layer1.Raster as IRasterProps;
                IRasterProps raster2005 = layer2.Raster as IRasterProps;

                double width = raster1996.Width;
                double height = raster1996.Height;


                //根据影像设定分辨率
                double cellWidth = raster1996.MeanCellSize().X;
                double cellHeight = raster1996.MeanCellSize().Y;

                IPnt pntSize = new PntClass();
                pntSize.SetCoords(cellWidth, cellHeight);
                IPixelBlock pixelBlock = layer1.Raster.CreatePixelBlock(pntSize);

                IRaster2 rst2First = layer1.Raster as IRaster2;
                IRaster2 rst2Sec = layer2.Raster as IRaster2;

                int count = 0;
                string str = "影像分辨率为:" + height + "×" + width;
                for (int i = 0; i < height; ++i)
                {
                    //DateTime dt1 = DateTime.Now;
                    for (int j = 0; j < width; ++j)
                    {
                        count++;
                        if (count % 50 == 0)
                        {
                            this.statusBarXY.Text = str + ";正在处理第" + i + "行,第" + j + "列";
                            System.Windows.Forms.Application.DoEvents();
                        }

                        int value1 = Function.GetPixelValue(j, i, layer1, pixelBlock);
                        if (value1 <= 0)
                            continue;
                        int value2 = Function.GetPixelValue(j, i, layer2, pixelBlock);
                        if (value2 <= 0)
                            continue;
                        if (landCount.ContainsKey(value1))
                        {
                            landCount[value1]++;
                        }
                        else
                        {
                            landCount.Add(value1, 1);
                        }
                        if (value1 == value2)
                        {
                            if (RightCount.ContainsKey(value1))
                            {
                                RightCount[value1]++;
                            }
                            else

                                RightCount.Add(value1, 1);
                        }
                    }
                }
            }

            StreamWriter sw = new StreamWriter(Application.StartupPath + "\\精度.txt");
            foreach (var key in landCount.Keys)
            {
                double a = (double)RightCount[key] / (double)landCount[key];
                string str = key + ":" + a;
                sw.WriteLine(str);
            }
            sw.Close();

            MessageBox.Show("计算结束。");
        }


        private void 道格拉斯普客法线ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            for (int k = 0; k < axMapControl1.LayerCount; k++)
            {

                IFeatureLayer featureLayer = axMapControl1.get_Layer(k) as IFeatureLayer;
                int FeatureCount = featureLayer.FeatureClass.FeatureCount(null); ;
                for (int i = 0; i < FeatureCount; i++)
                {
                    IFeature pFeature = featureLayer.FeatureClass.GetFeature(i);

                    //IPointCollection pCoor = pFeature.Shape as IPointCollection;
                    IPolyline polyline = pFeature.Shape as IPolyline;
                    polyline.Generalize(1);


                    pFeature.Store();

                    statusBarXY.Text = "共有：" + FeatureCount + "个要素，正在处理第" + i + "个";
                    Application.DoEvents();
                }

            }
            MessageBox.Show("处理完成！");
        }

        /****点到直线的距离***
  * 过点（x1,y1）和点（x2,y2）的直线方程为：KX -Y + (x2y1 - x1y2)/(x2-x1) = 0
  * 设直线斜率为K = (y2-y1)/(x2-x1),C=(x2y1 - x1y2)/(x2-x1)
  * 点P(x0,y0)到直线AX + BY +C =0DE 距离为：d=|Ax0 + By0 + C|/sqrt(A*A + B*B)
  * 点（x3,y3）到经过点（x1,y1）和点（x2,y2）的直线的最短距离为：
  * distance = |K*x3 - y3 + C|/sqrt(K*K + 1)
  */
        public double GetMinDistance(IPoint pt1, IPoint pt2, IPoint pt3)
        {
            double dis = 0;
            if (pt1.X == pt2.X)
            {
                dis = Math.Abs(pt3.X - pt1.X);
                return dis;
            }
            double lineK = (pt2.Y - pt1.Y) / (pt2.X - pt1.X);
            double lineC = (pt2.X * pt1.Y - pt1.X * pt2.Y) / (pt2.X - pt1.X);
            dis = Math.Abs(lineK * pt3.X - pt3.Y + lineC) / (Math.Sqrt(lineK * lineK + 1));
            return dis;
        }

        /// <summary>
        /// 根据最小值，平均值，标准差，最大值确定警限
        /// </summary>
        /// <returns></returns>
        public double[] GetWarnLevelValue(double min, double avg, double delt, double max)
        {
            double[] result = new double[6];
            //result[0] = min;
            result[1] = avg;
            result[2] = avg + 0.5 * delt;
            result[3] = avg + 1.5 * delt;
            result[4] = avg + 2.5 * delt;
            //result[2] = avg + delt;
            //result[3] = avg + 2 * delt;
            //result[4] = avg + 3 * delt;
            if ((avg + 4 * delt) > max)
            {
                result[5] = avg + 4 * delt;
            }
            else
                result[5] = max;
            return result;
        }

        public double[] GetWarnLevelValue(double min, double mid, double max)
        {
            double[] result = new double[6];
            result[0] = max;
            result[1] = mid;
            double length = (mid - min) / 4;
            result[2] = mid - length;
            result[3] = mid - 2 * length;
            result[4] = mid - 3 * length;
            result[5] = min;

            return result;
        }
        public WarnLevel[,] ReadWarnLevel(string filepath, ref double[] v)
        {
            int rows = 20;
            WarnLevel[,] result = new WarnLevel[rows, 5];
            return result;
        }
        /// <summary>
        /// 输出预警结果
        /// </summary>
        /// <param name="a"></param>
        /// <returns></returns>
        public string GetResult(double[] a)
        {
            int index = GetMaxIndex(a);
            if (index == 0)
                return "无警";
            else if (index == 1)
                return "轻警";
            else if (index == 2)
                return "中警";
            else if (index == 3)
                return "重警";
            else
                return "巨警";

        }

        public string GetResult(double a)
        {

            if (a >= 80)
                return "无警";
            else if (a >= 60)
                return "轻警";
            else if (a >= 40)
                return "中警";
            else if (a >= 20)
                return "重警";
            else
                return "巨警";

        }
        public int GetMaxIndex(double[] a)
        {
            double max = a[0];
            int index = 0;
            for (int i = 1; i < a.Length; i++)
            {
                if (a[i] > max)
                {
                    max = a[i];
                    index = i;
                }

            }
            return index;
        }
        public int GetMaxIndex(double a)
        {
            if (a >= 80)
                return 1;
            else if (a >= 60)
                return 2;
            else if (a >= 40)
                return 3;
            else if (a >= 20)
                return 4;
            else
                return 5;


        }
        /// <summary>
        /// 求关联系数
        /// </summary>
        /// <param name="v">待评价区域量值</param>
        /// <param name="warnlevel">经典域矩阵</param>
        /// <param name="i">第i个指标</param>
        /// <param name="j">第j个评价等级</param>
        /// <returns></returns>
        public double CalCorrelation(double v, WarnLevel warnlevel, int i, int j)
        {
            double Correlation = 0;

            singleWarnLevel wlIJ = warnlevel.GetSingleLevel(j + 1);//warnlevel[i, j+1];
            singleWarnLevel wPI = GetPI(i, warnlevel);
            if (wlIJ.In(v))//v属于V0
            {
                Correlation = -P(v, wlIJ) / wlIJ.GetLength();
            }
            else//不属于
            {
                Correlation = P(v, wlIJ) / (P(v, wPI) - P(v, wlIJ));
            }
            return Correlation;
        }
        public double P(double v, singleWarnLevel wl)
        {
            double result = 0;
            result = Math.Abs(v - 0.5 * (wl.low + wl.high)) - 0.5 * (wl.high - wl.low);
            return result;
        }
        /// <summary>
        /// 求平均值
        /// </summary>
        /// <param name="Valist"></param>
        /// <returns></returns>
        public double average(double[] Valist)
        {
            double sum = 0;
            foreach (double d in Valist)
            {
                sum = sum + d;
            }
            double revl = sum / Valist.Length;
            return revl;
        }

        /// <summary>
        /// 求数据集合标准差
        /// </summary>
        /// <param name="ValList"></param>
        /// <returns></returns>
        public double stdev(double[] ValList)
        {
            double avg = average(ValList);
            double sumstdev = 0;
            foreach (double d in ValList)
            {
                sumstdev += (d - avg) * (d - avg);
            }
            sumstdev /= ValList.Length;
            double stdeval = System.Math.Sqrt(sumstdev);
            return stdeval;
        }

        public double GetScore(double value, WarnLevel wl)
        {
            double score = 0;
            for (int i = 0; i < 5; i++)
            {
                singleWarnLevel swl = wl.GetSingleLevel(i + 1);
                if (swl.In(value))
                {
                    double H = (5 - i) * 20;
                    double L = (5 - i - 1) * 20;
                    score = (H - L) / (swl.high - swl.low) * (value - swl.low) + L;
                    return score;
                }
            }

            return score;
        }
        /// <summary>
        /// 获取节域物元矩阵
        /// </summary>
        /// <param name="i"></param>
        /// <param name="warnLevel"></param>
        /// <returns></returns>
        public singleWarnLevel GetPI(int i, WarnLevel warnLevel)
        {

            return new singleWarnLevel(warnLevel.level0, warnLevel.level5);
        }

        //计算更指标值的权重
        public List<double> CalcWeight(List<RoadObject> roads, List<double> sum)
        {
            List<double> weight = new List<double>();
            int ObjectCount = roads.Count;
            int IndicatorCount = roads[0].Indicators.Count;
            List<List<double>> PIJ = new List<List<double>>();
            #region 无量纲化
            for (int i = 0; i < ObjectCount; i++)
            {
                List<Double> indicatorValue = new List<double>();
                for (int j = 0; j < IndicatorCount; j++)
                {
                    double x = roads[i].Indicators[j] / sum[j];
                    indicatorValue.Add(x);
                }
                PIJ.Add(indicatorValue);
            }
            #endregion

            #region 计算熵值
            List<double> shang = new List<double>();//熵值
            for (int i = 0; i < IndicatorCount; i++)
            {
                double x = 0;
                for (int j = 0; j < ObjectCount; j++)
                {
                    if (PIJ[j][i] > 0)
                        x += PIJ[j][i] * Math.Log(Math.E, PIJ[j][i]);
                }
                x *= -1;
                shang.Add(x);
            }
            #endregion

            #region 计算权重值
            List<double> tmph = new List<double>();
            double hsum = 0;
            for (int i = 0; i < IndicatorCount; i++)
            {
                double tmp = shang[i] / Math.Log(Math.E, ObjectCount);
                tmph.Add(tmp);
                hsum += (1 - tmp);
            }
            for (int i = 0; i < IndicatorCount; i++)
            {
                double wei = (1 - tmph[i]) / hsum;
                weight.Add(wei);
            }

            #endregion
            tmph.Clear();
            shang.Clear();
            for (int i = 0; i < IndicatorCount; i++)
            {
                PIJ[i].Clear();
            }
            PIJ.Clear();
            return weight;
        }

        private void 批量计算混合利用指数ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            for (int k = 0; k < m_mapControl.LayerCount; k++)
            {
                IRasterLayer layer = m_mapControl.get_Layer(k) as IRasterLayer;
                IRasterProps pRP1 = layer.Raster as IRasterProps;

                //影像的宽度和高度
                double width = pRP1.Width;
                double height = pRP1.Height;


                //根据影像设定分辨率
                double cellWidth = pRP1.MeanCellSize().X;
                double cellHeight = pRP1.MeanCellSize().Y;
                IPnt pntSize = new PntClass();
                pntSize.SetCoords(cellWidth, cellHeight);
                IPixelBlock pixelBlock = layer.Raster.CreatePixelBlock(pntSize);
                double cellSize = cellWidth;
                IRaster2 rst2First = layer.Raster as IRaster2;
                IRaster raster = layer.Raster;
                IRasterProps raster1996 = raster as IRasterProps;

                //   指定像素块大小
                IPnt pntBlock = new PntClass();
                pntBlock.X = raster1996.Width;
                pntBlock.Y = raster1996.Height;
                //创建一个光标以给定像素块大小
                IRasterCursor rasterCursor = rst2First.CreateCursorEx(pntBlock);
                //创建像素块
                IPixelBlock3 pixelBlock3 = null;
                int blockWidth = 0;
                int blockHeight = 0;
                pixelBlock3 = rasterCursor.PixelBlock as IPixelBlock3;
                blockWidth = pixelBlock3.Width;
                blockHeight = pixelBlock3.Height;

                //指定平面的像素的数组
                System.Array pixels = (System.Array)pixelBlock3.get_PixelData(0);
                System.Array tmppixels = (System.Array)pixels.Clone();
                double nodata = Convert.ToDouble(tmppixels.GetValue(0, 0));
                int times = 0;
                string str = "共有" + m_mapControl.LayerCount + "个图层，正在处理第" + (k + 1) + "个。";
                str += "当前影像分辨率为:" + height + "×" + width;
                for (int i = 0; i < height; ++i)
                {
                    for (int j = 0; j < width; ++j)
                    {
                        times++;
                        if (times >= 500)
                        {
                            //str += ";正在处理第" + i + "行,第" + j + "列";
                            this.statusBarXY.Text = str + ";正在处理第" + (i + 1) + "行,第" + (j + 1) + "列";
                            System.Windows.Forms.Application.DoEvents();
                            times = 0;
                        }

                        double value1 = Convert.ToDouble(tmppixels.GetValue(j, i));
                        if (value1 == nodata)
                            continue;
                        double result = AEFunction.ArcEngineFunction.Calc101X101(j, i, layer, pixelBlock, tmppixels);
                        pixels.SetValue((float)result, j, i);

                    }
                }

                IRasterEdit rasterEdit = rst2First as IRasterEdit;
                if (rasterEdit.CanEdit())
                {

                    pixelBlock3.set_PixelData(0, pixels);
                    IPnt pnt = null;
                    pnt = rasterCursor.TopLeft;
                    rasterEdit.Write(pnt, (IPixelBlock)pixelBlock3);
                    rasterEdit.Refresh();
                }

                System.Runtime.InteropServices.Marshal.ReleaseComObject(rasterEdit);
            }
        }

        private void toolStripMenuItem2_Click(object sender, EventArgs e)
        {
            int count = m_mapControl.LayerCount;
            for (int k = 0; k < count; k++)
            {

                string str = "共有:" + count + "个图层;正在处理第" + k + "个";
                this.statusBarXY.Text = str;
                IRasterLayer layer = m_mapControl.get_Layer(k) as IRasterLayer;
                Geoprocessor GP = new Geoprocessor();
                ESRI.ArcGIS.DataManagementTools.CopyRaster cr = new ESRI.ArcGIS.DataManagementTools.CopyRaster();
                cr.in_raster = layer;
                cr.out_rasterdataset = @"F:\s私活\预警数据\09-20\double.mdb\" + k.ToString();
                cr.pixel_type = "64_BIT";
                GP.Execute(cr, null);
            }
            MessageBox.Show("导出完成！");
        }

        private void clipByFeatureToRasterToolStripMenuItem_Click(object sender, EventArgs e)
        {
            List<ILayer> layer = new List<ILayer>();
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                layer.Add(axMapControl1.get_Layer(i));
            }
            AEFunction.UI.BatchClipByFeatureToRaster batchClip = new AEFunction.UI.BatchClipByFeatureToRaster(toolStripProgressBar1, axMapControl1, layer.ToArray());
            batchClip.ShowDialog();
        }

        private void 模糊可拓预警ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            OpenFileDialog ofd = new OpenFileDialog();
            //ofd.InitialDirectory = Application.StartupPath;
            ofd.Filter = "指标文件(*.xlsx)|*.xlsx|指标文件(*.xls)|*.xls|所有文件(*.*)|*.*";
            //SELECT DISTINCTROW Sum([Sheet1].[2017年扩展坡度构成指数标准化值]) FROM Sheet1;

            FolderBrowserDialog fbd = new FolderBrowserDialog();
            if (fbd.ShowDialog() != DialogResult.OK)
                return;
            string savePath = fbd.SelectedPath;
            if (ofd.ShowDialog() == DialogResult.OK)
            {
                object missing = System.Reflection.Missing.Value;
                string excelFileName = ofd.FileName;
                string strCon = " Provider =Microsoft.ACE.OleDb.12.0 ; Data Source =" + excelFileName + ";Extended Properties='Excel 12.0; HDR=Yes; IMEX=1'";
                OleDbConnection myConn = new OleDbConnection(strCon);
                string strCom = "select * from [Sheet1$]";
                OleDbDataAdapter myCommand = new OleDbDataAdapter(strCom, myConn);
                myConn.Open();
                List<RoadObject> AllRoad = new List<RoadObject>();
                List<double> sums = new List<double>();//记录每个指标所有值的总和值，顺序按EXCEL表来确定
                //List<double> avgs = new List<double>();//记录每个指标所有值的平均值，顺序按EXCEL表来确定
                //List<double> deltas = new List<double>();//记录每个指标的标准差，顺序按EXCEL表来确定
                List<double> mins = new List<double>();//记录每个指标的最小值，顺序按EXCEL表来确定
                List<double> maxs = new List<double>();//记录每个指标的最小值，顺序按EXCEL表来确定
                List<double> mids = new List<double>();//记录每个指标的2/3位置的值，顺序按EXCEL表来确定
                List<double> weight = new List<double>();//记录每个指标的权重

                OpenFileDialog weightDlg = new OpenFileDialog();
                weightDlg.Filter = "指标权重文件(*.txt)|*.txt|所有文件(*.*)|*.*";
                if (weightDlg.ShowDialog() == DialogResult.OK)
                {
                    StreamReader sr = new StreamReader(weightDlg.FileName);
                    string line;
                    while ((line = sr.ReadLine()) != null)
                    {
                        weight.Add(Convert.ToDouble(line.Split('_')[1]));
                    }
                }
                DataTable myDataTable = new DataTable();
                List<WarnLevel> warnLevels = new List<WarnLevel>();//记录每个指标的警限，顺序按EXCEL表来确定
                string indicatorName = "街道名称 ";
                #region 读取道路类，计算各个指标值的总和、标准差
                myCommand.Fill(myDataTable);
                for (int i = 1; i < myDataTable.Columns.Count; i++)
                {
                    string filedName = myDataTable.Columns[i].Caption;
                    indicatorName += filedName + " ";
                    //double ss1 = (double)myDataTable.Compute("Avg([" + filedName + "])", "");
                    //avgs.Add(ss1);
                    //double ss2 = (double)myDataTable.Compute("Sum([" + filedName + "])", "");
                    //sums.Add(ss2);
                    double ss3 = (double)myDataTable.Compute("max([" + filedName + "])", "");
                    maxs.Add(ss3);
                    double ss4 = (double)myDataTable.Compute("min([" + filedName + "])", "");
                    mins.Add(ss4);
                }
                for (int i = 0; i < myDataTable.Rows.Count; i++)
                {
                    string roadName = myDataTable.Rows[i][0] as string;
                    RoadObject road = new RoadObject(roadName);
                    for (int j = 1; j < myDataTable.Columns.Count; j++)
                    {
                        double values = (double)myDataTable.Rows[i][j];
                        road.Indicators.Add(values);

                    }
                    AllRoad.Add(road);
                }
                int roadCount = AllRoad.Count;
                for (int i = 0; i < AllRoad[0].Indicators.Count; i++)
                {
                    List<double> aa = new List<double>();
                    for (int j = 0; j < AllRoad.Count; j++)
                    {
                        aa.Add(AllRoad[j].Indicators[i]);
                    }
                    aa.Sort();
                    mids.Add(aa[(int)roadCount * 1 / 3]);
                    //deltas.Add(stdev(aa.ToArray()));
                    aa.Clear();
                    #region 警限的确定
                    double[] tmpWL = GetWarnLevelValue(mins[i], mids[i], maxs[i]);
                    //double[] tmpWL = GetWarnLevelValue(mins[i], avgs[i], deltas[i], maxs[i]);
                    warnLevels.Add(new WarnLevel(tmpWL[0], tmpWL[1], tmpWL[2], tmpWL[3], tmpWL[4], tmpWL[5]));
                    #endregion
                }
                #endregion

                if (weight.Count == 0)
                    weight = CalcWeight(AllRoad, sums);

                int IndicatorCount = AllRoad[0].Indicators.Count;
                FileStream fs1 = new FileStream(savePath + "\\分项预警指标.txt", FileMode.OpenOrCreate);
                StreamWriter sw1 = new StreamWriter(fs1);
                FileStream fs2 = new FileStream(savePath + "\\综合预警指标.txt", FileMode.OpenOrCreate);
                StreamWriter sw2 = new StreamWriter(fs2);
                sw1.WriteLine(indicatorName);
                foreach (RoadObject ro in AllRoad)
                {
                    sw1.Write(ro.Name + " ");
                    sw2.Write(ro.Name + " ");
                    //List<double[]> CorrelationArray = new List<double[]>();//关联度矩阵
                    List<double> IndicatorScores = new List<double>();//关联度矩阵
                    for (int i = 0; i < IndicatorCount; i++)
                    {
                        //double[] correlas = new double[5];
                        //for (int j = 0; j < 5; j++)
                        //{
                        //    correlas[j] = CalCorrelation(ro.Indicators[i], warnLevels[i], i, j);
                        //}
                        //CorrelationArray.Add(correlas);
                        //sw1.Write(GetResult(CorrelationArray[i]) + " ");

                        double score = 0; ;
                        //for (int j = 0; j < 5; j++)
                        //{
                        //    correlas[j] = CalCorrelation(ro.Indicators[i], warnLevels[i], i, j);
                        //}
                        score = GetScore(ro.Indicators[i], warnLevels[i]);
                        IndicatorScores.Add(score);
                        sw1.Write(GetResult(score) + " ");

                    }
                    sw1.WriteLine();
                    //double[] ComprehensiveCorrelation = new double[5];//综合关联度
                    //for (int i = 0; i < 5; i++)
                    //{
                    //    ComprehensiveCorrelation[i] = 0;
                    //    for (int j = 0; j < IndicatorCount; j++)
                    //    {
                    //        ComprehensiveCorrelation[i] += weight[j] * IndicatorScores[j];
                    //    }
                    //}
                    double totleScore = 0;
                    for (int j = 0; j < IndicatorCount; j++)
                    {
                        totleScore += weight[j] * IndicatorScores[j];
                    }
                    sw2.Write(GetResult(totleScore));
                    sw2.WriteLine();
                }
                sw2.Close();
                sw1.Close();
                fs1.Close();
                fs2.Close();
                myConn.Close();
                myConn.Dispose();
                myCommand.Dispose();
                MessageBox.Show("预警完成。");
            }
        }
        private void 模糊可拓预警栅格尺度ToolStripMenuItem_Click(object sender, EventArgs e)
        {


            FolderBrowserDialog fbd = new FolderBrowserDialog();
            if (fbd.ShowDialog() != DialogResult.OK)
                return;
            string savePath = fbd.SelectedPath;
            Geoprocessor tGp = new Geoprocessor();
            int indicatorCount = axMapControl1.LayerCount - 1;
            IRaster[] pRasters = new IRaster[indicatorCount];
            IRaster2[] pRaster2s = new IRaster2[indicatorCount];
            IRasterProps[] rasterProps = new IRasterProps[indicatorCount];
            Array[] pSafeArray = new Array[indicatorCount + 1];
            ILayer pointLayer = axMapControl1.get_Layer(0);

            List<RoadObject> AllRoad = new List<RoadObject>();//记录每个栅格
            List<WarnLevel> warnLevels = new List<WarnLevel>();//记录每个指标的警限，顺序按EXCEL表来确定
            List<double>[] values = new List<double>[indicatorCount];

            double[] tmpMax = new double[indicatorCount];
            double[] tmpMin = new double[indicatorCount];
            //double[] tmpSum = new double[indicatorCount];
            //double[] tmpAvg = new double[indicatorCount];
            //double[] tmpDev = new double[indicatorCount];
            double[] mids = new double[indicatorCount];//记录每个指标的2/3位置的值，顺序按EXCEL表来确定

            for (int i = 1; i < axMapControl1.LayerCount; i++)
            {
                pRasters[i - 1] = (axMapControl1.get_Layer(i) as IRasterLayer).Raster;
                pRaster2s[i - 1] = pRasters[i - 1] as IRaster2;
                rasterProps[i - 1] = (IRasterProps)pRasters[i - 1];

                IPnt pBlockSize = new Pnt();
                pBlockSize.SetCoords(rasterProps[i - 1].Width, rasterProps[i - 1].Height);
                //选取整个范围  
                IPixelBlock pPixelBlock = pRasters[i - 1].CreatePixelBlock(pBlockSize);
                //左上点坐标  
                IPnt tlp = new Pnt();
                tlp.SetCoords(0, 0);
                //读入栅格  
                IRasterBandCollection pRasterBands = pRasters[i - 1] as IRasterBandCollection;
                IRasterBand pRasterBand = pRasterBands.Item(0);
                IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
                pRawRixels.Read(tlp, pPixelBlock);
                //将PixBlock的值组成数组  
                pSafeArray[i - 1] = pPixelBlock.get_SafeArray(0) as System.Array;
                object ndv = rasterProps[i - 1].NoDataValue;

                #region 读取平均值等
                //ESRI.ArcGIS.DataManagementTools.GetRasterProperties GpGetRasterProp = new ESRI.ArcGIS.DataManagementTools.GetRasterProperties(axMapControl1.get_Layer(i) as IRasterLayer);
                //GpGetRasterProp.property_type = "MEAN";
                //tGp.Execute(GpGetRasterProp, null);
                //tmpAvg[i-1]= (double)GpGetRasterProp.property;

                //GpGetRasterProp.property_type = "MINIMUM";
                //tGp.Execute(GpGetRasterProp, null);
                //tmpMin[i-1] = (double)GpGetRasterProp.property;

                //GpGetRasterProp.property_type = "MAXIMUM";
                //tGp.Execute(GpGetRasterProp, null);
                //tmpMax[i - 1] = (double)GpGetRasterProp.property;


                //GpGetRasterProp.property_type = "STD";
                //tGp.Execute(GpGetRasterProp, null);
                //tmpDev[i - 1] = (double)GpGetRasterProp.property;
                #endregion

            }
            pSafeArray[indicatorCount] = (Array)pSafeArray[0].Clone();
            int width = rasterProps[0].Width;
            int height = rasterProps[0].Height;
            //分项预警指标


            double cellsize = rasterProps[0].MeanCellSize().X;

            ISpatialReference sr = rasterProps[0].SpatialReference;
            IPoint originPoint = rasterProps[0].Extent.LowerLeft;

            //综合预警指标
            //Array pTotalPixels = Array.CreateInstance(typeof(int), width, height);

            int cellcounts = (pointLayer as IFeatureLayer).FeatureClass.FeatureCount(null);

            List<double> weight = new List<double>();//记录每个指标的权重

            OpenFileDialog weightDlg = new OpenFileDialog();
            weightDlg.Filter = "指标权重文件(*.txt)|*.txt|所有文件(*.*)|*.*";
            if (weightDlg.ShowDialog() == DialogResult.OK)
            {
                StreamReader sr1 = new StreamReader(weightDlg.FileName);
                string line;
                while ((line = sr1.ReadLine()) != null)
                {
                    weight.Add(Convert.ToDouble(line.Split('_')[1]));
                }
            }


            for (int i = 0; i < indicatorCount; i++)
            {
                values[i] = new List<double>();
            }
            for (int i = 0; i < cellcounts; i++)
            {
                //System.Windows.Forms.Application.DoEvents();
                IFeature f = (pointLayer as IFeatureLayer).FeatureClass.GetFeature(i);
                IPoint p = f.Shape as IPoint;
                int y = pRaster2s[0].ToPixelRow(p.Y);
                int x = pRaster2s[0].ToPixelColumn(p.X);
                RoadObject ro = new RoadObject(y + "_" + x);
                for (int j = 0; j < indicatorCount; j++)
                {
                    //int y = pRaster2s[j].ToPixelRow(p.Y);
                    //int x = pRaster2s[j].ToPixelColumn(p.X);
                    object oo = pSafeArray[j].GetValue(x, y);
                    double value = -1;
                    double.TryParse(oo.ToString(), out value);
                    ro.Indicators.Add(value);

                    values[j].Add(value);
                }

                AllRoad.Add(ro);
            }

            for (int i = 0; i < indicatorCount; i++)
            {
                values[i].Clear();
            }
            Array.Clear(values, 0, values.Length);
            Array.Clear(pSafeArray[0], 0, pSafeArray[0].Length);
            int roadCount = AllRoad.Count;
            //警限的确定
            for (int i = 0; i < AllRoad[0].Indicators.Count; i++)
            {
                values[i].Sort();
                tmpMin[i] = values[i][0];
                mids[i] = values[i][(int)roadCount * 1 / 3];
                tmpMax[i] = values[i][roadCount - 1];
                //deltas.Add(stdev(aa.ToArray()));
                values[i].Clear();
                double[] tmpWL = GetWarnLevelValue(tmpMin[i], mids[i], tmpMax[i]);
                warnLevels.Add(new WarnLevel(tmpWL[0], tmpWL[1], tmpWL[2], tmpWL[3], tmpWL[4], tmpWL[5]));
            }
            //List<double> sums = new List<double>(tmpSum);
            //if (weight.Count == 0)
            //    weight = CalcWeight(AllRoad, sums);
            foreach (RoadObject ro in AllRoad)
            {

                int row = Convert.ToInt32(ro.Name.Split('_')[0]);
                int col = Convert.ToInt32(ro.Name.Split('_')[1]);
                List<double> IndicatorScores = new List<double>();//关联度矩阵
                for (int i = 0; i < indicatorCount; i++)
                {
                    double score = 0; ;

                    score = GetScore(ro.Indicators[i], warnLevels[i]);
                    IndicatorScores.Add(score);

                    int result = GetMaxIndex(score);
                    pSafeArray[i].SetValue((float)result, col, row);
                    //sw1.Write(GetResult(CorrelationArray[i]) + " ");
                }
                //sw1.WriteLine();
                double totleScore = 0;//综合关联度
                for (int i = 0; i < 5; i++)
                {
                    totleScore += weight[i] * IndicatorScores[i];
                }
                int totalresult = GetMaxIndex(totleScore);
                pSafeArray[indicatorCount].SetValue((float)totalresult, col, row);
            }
            IRasterWorkspace2 rws = AEFunction.ArcEngineFunction.CreateRasterWorkspace(savePath);
            for (int i = 0; i < indicatorCount; i++)
            {
                string LayerName = axMapControl1.get_Layer(i + 1).Name;
                if (!LayerName.Contains(".tif"))
                {
                    LayerName += ".tif";
                }
                ArcEngineFunction.CreateFileRasterDataset(rws, "预警结果" + LayerName, originPoint, width, height, cellsize, rstPixelType.PT_FLOAT, sr, pSafeArray[i]);
            }
            ArcEngineFunction.CreateFileRasterDataset(rws, "综合预警指标.tif", originPoint, width, height, cellsize, rstPixelType.PT_FLOAT, sr, pSafeArray[indicatorCount]);


            System.Runtime.InteropServices.Marshal.ReleaseComObject(rws);
            MessageBox.Show("预警完成。");
        }
        private void ExportRasterToTxtToolStripMenuItem_Click(object sender, EventArgs e)
        {

            FolderBrowserDialog fbd = new FolderBrowserDialog();
            if (fbd.ShowDialog() != DialogResult.OK)
                return;
            string savePath = fbd.SelectedPath;
            int indicatorCount = axMapControl1.LayerCount - 1;
            ILayer pointLayer = axMapControl1.get_Layer(0);
            int cellcounts = (pointLayer as IFeatureLayer).FeatureClass.FeatureCount(null);
            List<RoadObject> AllRoad = new List<RoadObject>();//记录每个栅格
            List<WarnLevel> warnLevels = new List<WarnLevel>();//记录每个指标的警限，顺序按EXCEL表来确定
            List<double>[] values = new List<double>[indicatorCount];
            for (int i = 0; i < indicatorCount; i++)
            {
                values[i] = new List<double>();
            }
            for (int i = 1; i < axMapControl1.LayerCount; i++)
            {
                ILayer layer = axMapControl1.get_Layer(i);
                StreamWriter sw = new StreamWriter(savePath + "//" + i + layer.Name.Split('.')[0] + ".txt");
                string stastuStr = "共" + (axMapControl1.LayerCount - 1) + "个图层，正在读取第" + i + "个。  ";
                IRaster pRasters = (layer as IRasterLayer).Raster;
                IRaster2 pRaster2s = pRasters as IRaster2;
                IRasterProps rasterProps = (IRasterProps)pRasters;

                IPnt pBlockSize = new Pnt();
                pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);
                //选取整个范围  
                IPixelBlock pPixelBlock = pRasters.CreatePixelBlock(pBlockSize);
                //左上点坐标  
                IPnt tlp = new Pnt();
                tlp.SetCoords(0, 0);
                //读入栅格  
                IRasterBandCollection pRasterBands = pRasters as IRasterBandCollection;
                IRasterBand pRasterBand = pRasterBands.Item(0);
                IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
                pRawRixels.Read(tlp, pPixelBlock);
                //将PixBlock的值组成数组  
                Array pSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
                // object ndv = rasterProps[i - 1].NoDataValue;

                int times = 0;
                for (int j = 0; j < cellcounts; j++)
                {
                    times++;
                    if (times >= 2000)
                    {
                        System.Windows.Forms.Application.DoEvents();
                        statusBarXY.Text = stastuStr + (j + 1) + "/" + cellcounts;
                        times = 0;
                    }
                    IFeature f = (pointLayer as IFeatureLayer).FeatureClass.GetFeature(j);
                    IPoint p = f.Shape as IPoint;
                    int y = pRaster2s.ToPixelRow(p.Y);
                    int x = pRaster2s.ToPixelColumn(p.X);
                    double value = Convert.ToDouble(pSafeArray.GetValue(x, y));
                    sw.WriteLine(value);
                    //values[i-1].Add(value);
                    //if (i == 1)
                    //{
                    //    RoadObject ro = new RoadObject(y + "_" + x);
                    //    ro.Indicators.Add(value);
                    //    AllRoad.Add(ro);
                    //}
                    //else
                    //{
                    //    AllRoad[j].Indicators.Add(value);
                    //}
                }
                sw.Close();
                GC.Collect();


            }
        }
        private void 人均绿地面积ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            string name = axMapControl1.get_Layer(0).Name;
            Regex re = new Regex(@"\d.");
            string year = re.Match(name).Value;
            ILayer pointLayer = ArcEngineFunction.GetLayerByName(axMapControl1, year + "points");
            IRaster ldRasterLayer12 = (ArcEngineFunction.GetLayerByName(axMapControl1, year + "绿地面积.tif") as IRasterLayer).Raster;//绿地对应的图层
            IRaster pRJRasters = (ArcEngineFunction.GetLayerByName(axMapControl1, year + "人均绿地面积.tif") as IRasterLayer).Raster;
            IRaster pPOPRasters = (ArcEngineFunction.GetLayerByName(axMapControl1, year + "pop.tif") as IRasterLayer).Raster;//人口对应的点图层
            int cellcounts = (pointLayer as IFeatureLayer).FeatureClass.FeatureCount(null);
            //IFeatureSelection pLDFeatureSelection = ldPointLayer as IFeatureSelection;
            //IFeatureSelection pPOPFeatureSelection = popPointLayer as IFeatureSelection;
            #region  人均绿地面积
            IRaster2 pRaster2s = pRJRasters as IRaster2;
            IRasterProps rasterProps = (IRasterProps)pRJRasters;
            IPnt pBlockSize = new Pnt();
            pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock = pRJRasters.CreatePixelBlock(pBlockSize);
            //左上点坐标  
            IPnt tlp = new Pnt();
            tlp.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands = pRJRasters as IRasterBandCollection;
            IRasterBand pRasterBand = pRasterBands.Item(0);
            IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
            pRawRixels.Read(tlp, pPixelBlock);
            //将PixBlock的值组成数组  
            Array pRJSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
            rasterProps.NoDataValue = pRJSafeArray.GetValue(0, 0);
            #endregion
            #region  12年绿地面积
            IRaster2 pRaster2s12 = ldRasterLayer12 as IRaster2;
            IRasterProps rasterProps12 = (IRasterProps)ldRasterLayer12;
            IPnt pBlockSize12 = new Pnt();
            pBlockSize12.SetCoords(rasterProps12.Width, rasterProps12.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock12 = ldRasterLayer12.CreatePixelBlock(pBlockSize12);
            //左上点坐标  
            IPnt tlp12 = new Pnt();
            tlp12.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands12 = ldRasterLayer12 as IRasterBandCollection;
            IRasterBand pRasterBand12 = pRasterBands12.Item(0);
            IRawPixels pRawRixels12 = pRasterBands12.Item(0) as IRawPixels;
            pRawRixels12.Read(tlp12, pPixelBlock12);
            //将PixBlock的值组成数组  
            Array pldSafeArray12 = pPixelBlock12.get_SafeArray(0) as System.Array;
            int LDHeight12 = pldSafeArray12.GetLength(1);
            int LDWidth12 = pldSafeArray12.GetLength(0);
            #endregion
            #region 人口图层
            IRaster2 pRaster2s1 = pPOPRasters as IRaster2;
            IRasterProps rasterProps1 = (IRasterProps)pPOPRasters;
            IPnt pBlockSize1 = new Pnt();
            pBlockSize1.SetCoords(rasterProps1.Width, rasterProps1.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock1 = pPOPRasters.CreatePixelBlock(pBlockSize1);
            //左上点坐标  
            IPnt tlp1 = new Pnt();
            tlp1.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands1 = pPOPRasters as IRasterBandCollection;
            IRasterBand pRasterBand1 = pRasterBands1.Item(0);
            IRawPixels pRawRixels1 = pRasterBands1.Item(0) as IRawPixels;
            pRawRixels1.Read(tlp1, pPixelBlock1);
            //将PixBlock的值组成数组  
            Array pPOPSafeArray1 = pPixelBlock1.get_SafeArray(0) as System.Array;
            int popHeight = pPOPSafeArray1.GetLength(1);
            int popWidth = pPOPSafeArray1.GetLength(0);
            #endregion
            int times = 0;
            string statusBarStr = "";
            for (int i = 0; i < cellcounts; i++)
            {
                times++;
                if (times >= 500)
                {
                    Application.DoEvents();
                    statusBarStr = i + "//" + cellcounts;
                    this.statusBarXY.Text = statusBarStr;
                    times = 0;
                }
                IFeature f = (pointLayer as IFeatureLayer).FeatureClass.GetFeature(i);
                IPoint p = f.Shape as IPoint;
                double popSum = 0;
                double LDArea = 0;// (double)pLDFeatureSelection.SelectionSet.Count * 900;

                int tmpy = pRaster2s12.ToPixelRow(p.Y);
                int tmpx = pRaster2s12.ToPixelColumn(p.X);
                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= LDWidth12 || k >= LDHeight12)
                            continue;
                        double value = 0;
                        object oo = pldSafeArray12.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value == 0.0 || value == 255)
                            continue;
                        LDArea++;
                    }
                }
                LDArea *= 900;
                tmpy = pRaster2s1.ToPixelRow(p.Y);
                tmpx = pRaster2s1.ToPixelColumn(p.X);

                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= popWidth || k >= popHeight)
                            continue;
                        double value = 0;
                        object oo = pPOPSafeArray1.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value < -200)
                            continue;
                        popSum += value;
                    }
                }
                popSum *= 0.0009;
                if (popSum < 1)
                    popSum = 1;
                double result = LDArea / popSum;
                if (result > 80.953)
                    result = 80.953;
                int y = pRaster2s.ToPixelRow(p.Y);
                int x = pRaster2s.ToPixelColumn(p.X);

                pRJSafeArray.SetValue((float)result, x, y);

            }
            pPixelBlock.set_SafeArray(0, pRJSafeArray);
            IRasterEdit rasterEdit = pRJRasters as IRasterEdit;
            rasterEdit.Write(tlp, pPixelBlock);
            rasterEdit.Refresh();
            MessageBox.Show("处理完成。");
        }

        private void 人均绿地面积降速ToolStripMenuItem_Click(object sender, EventArgs e)
        {


            ILayer pointLayer = ArcEngineFunction.GetLayerByName(axMapControl1, "12points");
            IRaster ldRasterLayer12 = (ArcEngineFunction.GetLayerByName(axMapControl1, "12绿地面积.tif") as IRasterLayer).Raster;//12绿地对应的图层
            IRaster ldRasterLayer09 = (ArcEngineFunction.GetLayerByName(axMapControl1, "09绿地面积.tif") as IRasterLayer).Raster;//09绿地对应的图层
            //ILayer popPointLayer = ArcEngineFunction.GetLayerByName(axMapControl1, "12pop_points");
            IRaster pRJRasters = (ArcEngineFunction.GetLayerByName(axMapControl1, "12人均绿地面积降速.tif") as IRasterLayer).Raster;
            IRaster pPOPRasters = (ArcEngineFunction.GetLayerByName(axMapControl1, "SPOP2012_Resample1.tif") as IRasterLayer).Raster;//人口对应的点图层
            IRaster pPOPRasters2 = (ArcEngineFunction.GetLayerByName(axMapControl1, "SPOP2009_Resample1.tif") as IRasterLayer).Raster;//人口对应的点图层
            int cellcounts = (pointLayer as IFeatureLayer).FeatureClass.FeatureCount(null);

            #region 人均绿地面积降速
            IRaster2 pRaster2s = pRJRasters as IRaster2;
            IRasterProps rasterProps = (IRasterProps)pRJRasters;
            IPnt pBlockSize = new Pnt();
            pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);

            double dX = rasterProps.MeanCellSize().X; //栅格的宽度
            double dY = rasterProps.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock = pRJRasters.CreatePixelBlock(pBlockSize);
            //左上点坐标  
            IPnt tlp = new Pnt();
            tlp.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands = pRJRasters as IRasterBandCollection;
            IRasterBand pRasterBand = pRasterBands.Item(0);
            IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
            pRawRixels.Read(tlp, pPixelBlock);
            //将PixBlock的值组成数组  
            Array pRJSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
            #endregion
            #region  12年绿地面积
            IRaster2 pRaster2s12 = ldRasterLayer12 as IRaster2;
            IRasterProps rasterProps12 = (IRasterProps)ldRasterLayer12;
            IPnt pBlockSize12 = new Pnt();
            pBlockSize12.SetCoords(rasterProps12.Width, rasterProps12.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock12 = ldRasterLayer12.CreatePixelBlock(pBlockSize12);
            //左上点坐标  
            IPnt tlp12 = new Pnt();
            tlp12.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands12 = ldRasterLayer12 as IRasterBandCollection;
            IRasterBand pRasterBand12 = pRasterBands12.Item(0);
            IRawPixels pRawRixels12 = pRasterBands12.Item(0) as IRawPixels;
            pRawRixels12.Read(tlp12, pPixelBlock12);
            //将PixBlock的值组成数组  
            Array pldSafeArray12 = pPixelBlock12.get_SafeArray(0) as System.Array;
            int LDHeight12 = pldSafeArray12.GetLength(1);
            int LDWidth12 = pldSafeArray12.GetLength(0);
            #endregion

            #region  09年绿地面积
            IRaster2 pRaster2s09 = ldRasterLayer09 as IRaster2;
            IRasterProps rasterProps09 = (IRasterProps)ldRasterLayer09;
            IPnt pBlockSize09 = new Pnt();
            pBlockSize09.SetCoords(rasterProps09.Width, rasterProps09.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock09 = ldRasterLayer09.CreatePixelBlock(pBlockSize09);
            //左上点坐标  
            IPnt tlp09 = new Pnt();
            tlp09.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands09 = ldRasterLayer09 as IRasterBandCollection;
            IRasterBand pRasterBand09 = pRasterBands09.Item(0);
            IRawPixels pRawRixels09 = pRasterBands09.Item(0) as IRawPixels;
            pRawRixels09.Read(tlp09, pPixelBlock09);
            //将PixBlock的值组成数组  
            Array pldSafeArray09 = pPixelBlock09.get_SafeArray(0) as System.Array;
            int LDHeight09 = pldSafeArray09.GetLength(1);
            int LDWidth09 = pldSafeArray09.GetLength(0);
            #endregion
            #region 12人口图层
            IRaster2 pRaster2s1 = pPOPRasters as IRaster2;
            IRasterProps rasterProps1 = (IRasterProps)pPOPRasters;
            IPnt pBlockSize1 = new Pnt();
            pBlockSize1.SetCoords(rasterProps1.Width, rasterProps1.Height);


            double dX1 = rasterProps1.MeanCellSize().X; //栅格的宽度
            double dY1 = rasterProps1.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock1 = pPOPRasters.CreatePixelBlock(pBlockSize1);
            //左上点坐标  
            IPnt tlp1 = new Pnt();
            tlp1.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands1 = pPOPRasters as IRasterBandCollection;
            IRasterBand pRasterBand1 = pRasterBands1.Item(0);
            IRawPixels pRawRixels1 = pRasterBands1.Item(0) as IRawPixels;
            pRawRixels1.Read(tlp1, pPixelBlock1);
            //将PixBlock的值组成数组  
            Array pPOPSafeArray1 = pPixelBlock1.get_SafeArray(0) as System.Array;
            int popHeight = pPOPSafeArray1.GetLength(1);
            int popWidth = pPOPSafeArray1.GetLength(0);
            #endregion

            #region 09人口图层
            IRaster2 pRaster2s2 = pPOPRasters2 as IRaster2;
            IRasterProps rasterProps2 = (IRasterProps)pPOPRasters2;
            IPnt pBlockSize2 = new Pnt();
            pBlockSize2.SetCoords(rasterProps2.Width, rasterProps2.Height);


            double dX2 = rasterProps2.MeanCellSize().X; //栅格的宽度
            double dY2 = rasterProps2.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock2 = pPOPRasters2.CreatePixelBlock(pBlockSize2);
            //左上点坐标  
            IPnt tlp2 = new Pnt();
            tlp2.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands2 = pPOPRasters2 as IRasterBandCollection;
            IRasterBand pRasterBand2 = pRasterBands2.Item(0);
            IRawPixels pRawRixels2 = pRasterBands2.Item(0) as IRawPixels;
            pRawRixels2.Read(tlp2, pPixelBlock2);
            //将PixBlock的值组成数组  
            Array pPOPSafeArray2 = pPixelBlock2.get_SafeArray(0) as System.Array;
            int popHeight2 = pPOPSafeArray2.GetLength(1);
            int popWidth2 = pPOPSafeArray2.GetLength(0);
            #endregion
            string statusBarStr = "";
            for (int i = 0; i < cellcounts; i++)
            {
                //Application.DoEvents();
                statusBarStr = i + "//" + cellcounts;
                this.statusBarXY.Text = statusBarStr;
                IFeature f = (pointLayer as IFeatureLayer).FeatureClass.GetFeature(i);
                IPoint p = f.Shape as IPoint;

                double popSum = 0;
                double popSum09 = 0;
                double LDArea12 = 0;
                double LDArea09 = 0;

                int tmpy = pRaster2s12.ToPixelRow(p.Y);
                int tmpx = pRaster2s12.ToPixelColumn(p.X);
                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= LDWidth12 || k >= LDHeight12)
                            continue;
                        double value = 0;
                        object oo = pldSafeArray12.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value == 0.0)
                            continue;
                        LDArea12++;
                    }
                }
                //LDArea *= 900;

                tmpy = pRaster2s09.ToPixelRow(p.Y);
                tmpx = pRaster2s09.ToPixelColumn(p.X);
                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= LDWidth09 || k >= LDHeight09)
                            continue;
                        double value = 0;
                        object oo = pldSafeArray09.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value == 0.0)
                            continue;
                        LDArea09++;
                    }
                }
                //LDArea *= 900;


                tmpy = pRaster2s1.ToPixelRow(p.Y);
                tmpx = pRaster2s1.ToPixelColumn(p.X);

                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= popWidth || k >= popHeight)
                            continue;
                        double value = 0;
                        object oo = pPOPSafeArray1.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value < -200)
                            continue;
                        popSum += value;
                    }
                }
                popSum *= 0.0009;
                if (popSum < 1)
                    popSum = 1;

                int y = pRaster2s.ToPixelRow(p.Y);
                int x = pRaster2s.ToPixelColumn(p.X);
                if (LDArea09 == 0)
                {
                    float tmp = (float)((-LDArea12 * 900) / popSum);
                    pRJSafeArray.SetValue(tmp, x, y);
                    continue;
                }
                tmpy = pRaster2s2.ToPixelRow(p.Y);
                tmpx = pRaster2s2.ToPixelColumn(p.X);

                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= popWidth2 || k >= popHeight2)
                            continue;
                        double value = 0;
                        object oo = pPOPSafeArray2.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value < -200)
                            continue;
                        popSum09 += value;
                    }
                }
                popSum09 *= 0.0009;
                if (popSum09 < 1)
                    popSum09 = 1;

                double result = (LDArea09 / popSum09 - LDArea12 / popSum) / (LDArea09 / popSum09);

                // double result = buildArea / landCount;
                pRJSafeArray.SetValue((float)result, x, y);
            }
            pPixelBlock.set_SafeArray(0, pRJSafeArray);
            IRasterEdit rasterEdit = pRJRasters as IRasterEdit;
            rasterEdit.Write(tlp, pPixelBlock);
            rasterEdit.Refresh();
            MessageBox.Show("处理完成。");
        }

        private void 建筑密度ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            ILayer pointLayer = ArcEngineFunction.GetLayerByName(axMapControl1, "12points");

            IRaster buildingLayer = (ArcEngineFunction.GetLayerByName(axMapControl1, "y2012.tif") as IRasterLayer).Raster;//建筑图层
            IRaster pRasters = (ArcEngineFunction.GetLayerByName(axMapControl1, "12建筑密度.tif") as IRasterLayer).Raster;
            int cellcounts = (pointLayer as IFeatureLayer).FeatureClass.FeatureCount(null);

            #region 12建筑密度
            IRaster2 pRaster2s = pRasters as IRaster2;
            IRasterProps rasterProps = (IRasterProps)pRasters;
            IPnt pBlockSize = new Pnt();
            pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);

            double dX = rasterProps.MeanCellSize().X; //栅格的宽度
            double dY = rasterProps.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock = pRasters.CreatePixelBlock(pBlockSize);
            //左上点坐标  
            IPnt tlp = new Pnt();
            tlp.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands = pRasters as IRasterBandCollection;
            IRasterBand pRasterBand = pRasterBands.Item(0);
            IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
            pRawRixels.Read(tlp, pPixelBlock);
            //将PixBlock的值组成数组  
            Array pSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;

            #endregion
            #region 建筑面积
            IRaster2 pRaster2s09 = buildingLayer as IRaster2;
            IRasterProps rasterProps09 = (IRasterProps)buildingLayer;
            IPnt pBlockSize09 = new Pnt();
            pBlockSize09.SetCoords(rasterProps09.Width, rasterProps09.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock09 = buildingLayer.CreatePixelBlock(pBlockSize09);
            //左上点坐标  
            IPnt tlp09 = new Pnt();
            tlp09.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands09 = buildingLayer as IRasterBandCollection;
            IRasterBand pRasterBand09 = pRasterBands09.Item(0);
            IRawPixels pRawRixels09 = pRasterBands09.Item(0) as IRawPixels;
            pRawRixels09.Read(tlp09, pPixelBlock09);
            //将PixBlock的值组成数组  
            Array pldSafeArray09 = pPixelBlock09.get_SafeArray(0) as System.Array;
            int LDHeight09 = pldSafeArray09.GetLength(1);
            int LDWidth09 = pldSafeArray09.GetLength(0);
            #endregion
            int height = pSafeArray.GetLength(1);
            int width = pSafeArray.GetLength(0);

            string statusBarStr = "";
            for (int i = 0; i < cellcounts; i++)
            {
                Application.DoEvents();
                statusBarStr = i + "//" + cellcounts;
                this.statusBarXY.Text = statusBarStr;
                IFeature f = (pointLayer as IFeatureLayer).FeatureClass.GetFeature(i);
                IPoint p = f.Shape as IPoint;


                double buildArea = 0;
                double landCount = 0;

                int y = pRaster2s.ToPixelRow(p.Y);
                int x = pRaster2s.ToPixelColumn(p.X);

                for (int k = y - 50; k <= y + 50; k++)
                {
                    for (int j = x - 50; j <= x + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= width || k >= height)
                            continue;
                        double value = Convert.ToDouble(pSafeArray.GetValue(j, k));
                        if (value == 0.0)
                            continue;
                        landCount++;
                    }
                }

                int tmpy = pRaster2s09.ToPixelRow(p.Y);
                int tmpx = pRaster2s09.ToPixelColumn(p.X);
                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= LDWidth09 || k >= LDHeight09)
                            continue;
                        double value = 0;
                        object oo = pldSafeArray09.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value <= 0.0)
                            continue;
                        buildArea++;
                    }
                }

                double result = buildArea / landCount;
                pSafeArray.SetValue((float)result, x, y);

            }

            pPixelBlock.set_SafeArray(0, pSafeArray);
            IRasterEdit rasterEdit = pRasters as IRasterEdit;
            rasterEdit.Write(tlp, pPixelBlock);
            rasterEdit.Refresh();
            MessageBox.Show("处理完成。");
        }

        private void 交通密度ToolStripMenuItem_Click(object sender, EventArgs e)
        {

            string name = axMapControl1.get_Layer(0).Name;
            Regex re = new Regex(@"\d.");
            string year = re.Match(name).Value;

            ILayer pointLayer = ArcEngineFunction.GetLayerByName(axMapControl1, year + "points");

            IRaster buildingLayer = (ArcEngineFunction.GetLayerByName(axMapControl1, year + "交通.tif") as IRasterLayer).Raster;//建筑图层
            IRaster pRasters = (ArcEngineFunction.GetLayerByName(axMapControl1, year + "交通密度.tif") as IRasterLayer).Raster;
            int cellcounts = (pointLayer as IFeatureLayer).FeatureClass.FeatureCount(null);

            #region 12建筑密度
            IRaster2 pRaster2s = pRasters as IRaster2;
            IRasterProps rasterProps = (IRasterProps)pRasters;
            IPnt pBlockSize = new Pnt();
            pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);

            double dX = rasterProps.MeanCellSize().X; //栅格的宽度
            double dY = rasterProps.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock = pRasters.CreatePixelBlock(pBlockSize);
            //左上点坐标  
            IPnt tlp = new Pnt();
            tlp.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands = pRasters as IRasterBandCollection;
            IRasterBand pRasterBand = pRasterBands.Item(0);
            IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
            pRawRixels.Read(tlp, pPixelBlock);
            //将PixBlock的值组成数组  
            Array pSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;

            #endregion
            #region 交通面积
            IRaster2 pRaster2s09 = buildingLayer as IRaster2;
            IRasterProps rasterProps09 = (IRasterProps)buildingLayer;
            IPnt pBlockSize09 = new Pnt();
            pBlockSize09.SetCoords(rasterProps09.Width, rasterProps09.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock09 = buildingLayer.CreatePixelBlock(pBlockSize09);
            //左上点坐标  
            IPnt tlp09 = new Pnt();
            tlp09.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands09 = buildingLayer as IRasterBandCollection;
            IRasterBand pRasterBand09 = pRasterBands09.Item(0);
            IRawPixels pRawRixels09 = pRasterBands09.Item(0) as IRawPixels;
            pRawRixels09.Read(tlp09, pPixelBlock09);
            //将PixBlock的值组成数组  
            Array pldSafeArray09 = pPixelBlock09.get_SafeArray(0) as System.Array;
            int LDHeight09 = pldSafeArray09.GetLength(1);
            int LDWidth09 = pldSafeArray09.GetLength(0);
            #endregion
            int height = pSafeArray.GetLength(1);
            int width = pSafeArray.GetLength(0);

            string statusBarStr = "";
            int times = 0;
            for (int i = 0; i < cellcounts; i++)
            {
                times++;
                if (times >= 500)
                {
                    Application.DoEvents();
                    statusBarStr = i + "//" + cellcounts;
                    this.statusBarXY.Text = statusBarStr;
                    times = 0;
                }

                IFeature f = (pointLayer as IFeatureLayer).FeatureClass.GetFeature(i);
                IPoint p = f.Shape as IPoint;


                double buildArea = 0;
                double landCount = 0;

                int y = pRaster2s.ToPixelRow(p.Y);
                int x = pRaster2s.ToPixelColumn(p.X);

                for (int k = y - 50; k <= y + 50; k++)
                {
                    for (int j = x - 50; j <= x + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= width || k >= height)
                            continue;
                        double value = Convert.ToDouble(pSafeArray.GetValue(j, k));
                        if (value == 0.0 || Math.Abs(value) > 20000000000)
                            continue;
                        landCount++;
                    }
                }

                int tmpy = pRaster2s09.ToPixelRow(p.Y);
                int tmpx = pRaster2s09.ToPixelColumn(p.X);
                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= LDWidth09 || k >= LDHeight09)
                            continue;
                        double value = 0;
                        object oo = pldSafeArray09.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value <= 0.0 || value == 255)
                            continue;
                        buildArea++;
                    }
                }

                double result = buildArea / landCount;
                pSafeArray.SetValue((float)result, x, y);

            }

            pPixelBlock.set_SafeArray(0, pSafeArray);
            IRasterEdit rasterEdit = pRasters as IRasterEdit;
            rasterEdit.Write(tlp, pPixelBlock);
            rasterEdit.Refresh();
            MessageBox.Show("处理完成。");

        }

        private void 修复人均绿地面积ToolStripMenuItem_Click(object sender, EventArgs e)
        {

            ILayer pointLayer = ArcEngineFunction.GetLayerByName(axMapControl1, "12points");
            IRaster pRJRasters = (ArcEngineFunction.GetLayerByName(axMapControl1, "12人均绿地面积.tif") as IRasterLayer).Raster;
            int cellcounts = (pointLayer as IFeatureLayer).FeatureClass.FeatureCount(null);
            //IFeatureSelection pLDFeatureSelection = ldPointLayer as IFeatureSelection;
            //IFeatureSelection pPOPFeatureSelection = popPointLayer as IFeatureSelection;
            #region  人均绿地面积
            IRaster2 pRaster2s = pRJRasters as IRaster2;
            IRasterProps rasterProps = (IRasterProps)pRJRasters;
            IPnt pBlockSize = new Pnt();
            pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock = pRJRasters.CreatePixelBlock(pBlockSize);
            //左上点坐标  
            IPnt tlp = new Pnt();
            tlp.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands = pRJRasters as IRasterBandCollection;
            IRasterBand pRasterBand = pRasterBands.Item(0);
            IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
            pRawRixels.Read(tlp, pPixelBlock);
            //将PixBlock的值组成数组  
            Array pRJSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
            #endregion

            string statusBarStr = "";
            for (int i = 0; i < cellcounts; i++)
            {
                Application.DoEvents();
                statusBarStr = i + "//" + cellcounts;
                this.statusBarXY.Text = statusBarStr;
                IFeature f = (pointLayer as IFeatureLayer).FeatureClass.GetFeature(i);
                IPoint p = f.Shape as IPoint;
                int y = pRaster2s.ToPixelRow(p.Y);
                int x = pRaster2s.ToPixelColumn(p.X);

                object testoo = pRJSafeArray.GetValue(x, y);
                double testresult = 0;
                double.TryParse(testoo.ToString(), out testresult);
                if (testresult < 80.953)
                    continue;
                if (testresult > 80.953)
                    pRJSafeArray.SetValue((float)80.953, x, y);

            }
            pPixelBlock.set_SafeArray(0, pRJSafeArray);
            IRasterEdit rasterEdit = pRJRasters as IRasterEdit;
            rasterEdit.Write(tlp, pPixelBlock);
            rasterEdit.Refresh();
            MessageBox.Show("处理完成。");
        }

        private void 正向指标标准化ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                IRasterLayer rasterLayer = axMapControl1.get_Layer(i) as IRasterLayer;
                string statusText = "共有：" + axMapControl1.LayerCount + "个图层，正在处理第" + (i + 1) + "个;";
                if (rasterLayer == null)
                    continue;
                Normalization(rasterLayer, 1, statusText);
            }
            MessageBox.Show("标准化完成。");
        }

        /// <summary>
        /// 1正向2负向3适度
        /// </summary>
        /// <param name="rasterLayer"></param>
        /// <param name="e"></param>
        /// <param name="statusText"></param>
        private void Normalization(IRasterLayer rasterLayer, int e, string statusText)
        {
            Geoprocessor tGp = new Geoprocessor();
            string LayerPath = rasterLayer.FilePath;
            statusBarXY.Text = statusText + "正在构建信息表。";
            ESRI.ArcGIS.DataManagementTools.CalculateStatistics rasterStatistics = new ESRI.ArcGIS.DataManagementTools.CalculateStatistics(LayerPath);
            tGp.Execute(rasterStatistics, null);


            IRaster pRasters = rasterLayer.Raster;
            IRaster2 pRaster2s = pRasters as IRaster2;
            IRasterProps rasterProps = (IRasterProps)pRasters;
            rasterProps.PixelType = rstPixelType.PT_FLOAT;
            IPnt pBlockSize = new Pnt();
            pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);

            double dX = rasterProps.MeanCellSize().X; //栅格的宽度
            double dY = rasterProps.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock = pRasters.CreatePixelBlock(pBlockSize);
            //左上点坐标  
            IPnt tlp = new Pnt();
            tlp.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands = pRasters as IRasterBandCollection;
            IRasterBand pRasterBand = pRasterBands.Item(0);
            IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
            pRawRixels.Read(tlp, pPixelBlock);
            //将PixBlock的值组成数组  
            Array pSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
            double NoDataValue = 0;
            //double.TryParse(rasterProps.NoDataValue.ToString(), out NoDataValue);
            double.TryParse(pSafeArray.GetValue(0, 0).ToString(), out NoDataValue);

            if (e == 1)//正向指标
            {
                double tmpMax = 0;
                double tmpMin = 0;

                ESRI.ArcGIS.DataManagementTools.GetRasterProperties GpGetRasterProp = new ESRI.ArcGIS.DataManagementTools.GetRasterProperties(rasterLayer);
                GpGetRasterProp.property_type = "MAXIMUM";
                tGp.Execute(GpGetRasterProp, null);
                tmpMax = (double)GpGetRasterProp.property;


                GpGetRasterProp.property_type = "MINIMUM";
                tGp.Execute(GpGetRasterProp, null);
                tmpMin = (double)GpGetRasterProp.property;
                double length = tmpMax - tmpMin;
                int times = 0;
                if (tmpMin * tmpMax < 0)
                {
                    //double a1 = tmpMax;
                    //double a2 = (tmpMax > (-tmpMin)) ? tmpMax : (-tmpMin);
                    #region 极差化
                    for (int k = 0; k < rasterProps.Height; k++)
                    {
                        for (int j = 0; j < rasterProps.Width; j++)
                        {
                            if (times++ >= 2000)
                            {
                                Application.DoEvents();
                                statusBarXY.Text = statusText + "影像分辨率为:" + rasterProps.Height + "×" + rasterProps.Width + ";正在处理第" + (k + 1) + "行,第" + (j + 1) + "列";
                                times = 0;
                            }
                            double value = 0;
                            double.TryParse(pSafeArray.GetValue(j, k).ToString(), out value);
                            if (value - NoDataValue < 0.0001)
                                continue;
                            if (value >= 0)
                                pSafeArray.SetValue((float)(value / tmpMax), j, k);
                            else
                                pSafeArray.SetValue((float)((value - tmpMin) * (-tmpMin) / (length * tmpMax)), j, k);
                        }
                    }
                    #endregion
                }
                else
                {
                    for (int k = 0; k < rasterProps.Height; k++)
                    {
                        for (int j = 0; j < rasterProps.Width; j++)
                        {
                            if (times++ >= 2000)
                            {
                                Application.DoEvents();
                                statusBarXY.Text = statusText + "影像分辨率为:" + rasterProps.Height + "×" + rasterProps.Width + ";正在处理第" + (k + 1) + "行,第" + (j + 1) + "列";
                                times = 0;
                            }
                            double value = 0;
                            double.TryParse(pSafeArray.GetValue(j, k).ToString(), out value);
                            if (value == NoDataValue)
                                continue;
                            float result = (float)((value - tmpMin) / length);
                            pSafeArray.SetValue(result, j, k);
                        }
                    }
                }

            }
            else if (e == 2)//负向指标
            {
                double tmpMax = 0;
                double tmpMin = 0;

                ESRI.ArcGIS.DataManagementTools.GetRasterProperties GpGetRasterProp = new ESRI.ArcGIS.DataManagementTools.GetRasterProperties(rasterLayer);
                GpGetRasterProp.property_type = "MAXIMUM";
                tGp.Execute(GpGetRasterProp, null);
                tmpMax = (double)GpGetRasterProp.property;
                GpGetRasterProp.property_type = "MINIMUM";
                tGp.Execute(GpGetRasterProp, null);
                tmpMin = (double)GpGetRasterProp.property;
                double length = tmpMax - tmpMin;
                int times = 0;
                if (tmpMin * tmpMax < 0)
                {
                    //double a1 = tmpMax;
                    //double a2 = (tmpMax > (-tmpMin)) ? tmpMax : (-tmpMin);
                    #region 极差化
                    for (int k = 0; k < rasterProps.Height; k++)
                    {
                        for (int j = 0; j < rasterProps.Width; j++)
                        {
                            if (times++ >= 2000)
                            {
                                Application.DoEvents();
                                statusBarXY.Text = statusText + "影像分辨率为:" + rasterProps.Height + "×" + rasterProps.Width + ";正在处理第" + (k + 1) + "行,第" + (j + 1) + "列";
                                times = 0;
                            }
                            double value = 0;
                            double.TryParse(pSafeArray.GetValue(j, k).ToString(), out value);
                            if (value - NoDataValue < 0.0001)
                                continue;
                            if (value >= 0)
                                pSafeArray.SetValue((float)(value / tmpMax), j, k);
                            else
                                pSafeArray.SetValue((float)((tmpMax - value) * (-tmpMin) / (length * tmpMax)), j, k);
                        }
                    }
                    #endregion
                }
                else
                {
                    for (int k = 0; k < rasterProps.Height; k++)
                    {
                        for (int j = 0; j < rasterProps.Width; j++)
                        {
                            if (times++ >= 2000)
                            {
                                Application.DoEvents();
                                statusBarXY.Text = statusText + "影像分辨率为:" + rasterProps.Height + "×" + rasterProps.Width + ";正在处理第" + (k + 1) + "行,第" + (j + 1) + "列";
                                times = 0;
                            }
                            double value = 0;
                            double.TryParse(pSafeArray.GetValue(j, k).ToString(), out value);
                            if (value - NoDataValue < 0.0001)
                                continue;
                            pSafeArray.SetValue((float)((tmpMax - value) / length), j, k);

                        }
                    }
                }

            }
            else if (e == 3)//适度指标
            {
                double mean = 0;
                ESRI.ArcGIS.DataManagementTools.GetRasterProperties GpGetRasterProp = new ESRI.ArcGIS.DataManagementTools.GetRasterProperties(rasterLayer);
                GpGetRasterProp.property_type = "MEAN";
                tGp.Execute(GpGetRasterProp, null);
                mean = (double)GpGetRasterProp.property;

                int times = 0;

                for (int k = 0; k < rasterProps.Height; k++)
                {
                    for (int j = 0; j < rasterProps.Width; j++)
                    {
                        if (times++ >= 2000)
                        {
                            Application.DoEvents();
                            statusBarXY.Text = statusText + "影像分辨率为:" + rasterProps.Height + "×" + rasterProps.Width + ";正在处理第" + (k + 1) + "行,第" + (j + 1) + "列";
                            times = 0;
                        }
                        double value = 0;
                        double.TryParse(pSafeArray.GetValue(j, k).ToString(), out value);
                        if (value - NoDataValue < 0.0001)
                            continue;
                        if (value >= mean)
                            pSafeArray.SetValue((float)(mean / value), j, k);
                        else
                            pSafeArray.SetValue((float)(value / mean), j, k);

                    }
                }

            }
            pPixelBlock.set_SafeArray(0, pSafeArray);
            IRasterEdit rasterEdit = pRasters as IRasterEdit;
            rasterEdit.Write(tlp, pPixelBlock);
            rasterEdit.Refresh();

            statusBarXY.Text = "正在重建信息表。";
            tGp.Execute(rasterStatistics, null);
        }

        private void 负向指标标准化ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                IRasterLayer rasterLayer = axMapControl1.get_Layer(i) as IRasterLayer;
                string statusText = "共有：" + axMapControl1.LayerCount + "个图层，正在处理第" + (i + 1) + "个;";
                if (rasterLayer == null)
                    continue;
                Normalization(rasterLayer, 2, statusText);
            }
            MessageBox.Show("标准化完成。");
        }

        private void 适度指标标准化ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                IRasterLayer rasterLayer = axMapControl1.get_Layer(i) as IRasterLayer;
                string statusText = "共有：" + axMapControl1.LayerCount + "个图层，正在处理第" + (i + 1) + "个;";
                if (rasterLayer == null)
                    continue;
                Normalization(rasterLayer, 3, statusText);
            }
            MessageBox.Show("标准化完成。");
        }

        private void 负向指标负值设为1ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                IRasterLayer rasterLayer = axMapControl1.get_Layer(i) as IRasterLayer;
                string statusText = "共有：" + axMapControl1.LayerCount + "个图层，正在处理第" + (i + 1) + "个;";
                if (rasterLayer == null)
                    continue;
                IRaster pRasters = rasterLayer.Raster;
                IRaster2 pRaster2s = pRasters as IRaster2;
                IRasterProps rasterProps = (IRasterProps)pRasters;

                IPnt pBlockSize = new Pnt();
                pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);

                double dX = rasterProps.MeanCellSize().X; //栅格的宽度
                double dY = rasterProps.MeanCellSize().Y; //栅格的高度
                //选取整个范围  
                IPixelBlock pPixelBlock = pRasters.CreatePixelBlock(pBlockSize);
                //左上点坐标  
                IPnt tlp = new Pnt();
                tlp.SetCoords(0, 0);
                //读入栅格  
                IRasterBandCollection pRasterBands = pRasters as IRasterBandCollection;
                IRasterBand pRasterBand = pRasterBands.Item(0);
                IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
                pRawRixels.Read(tlp, pPixelBlock);
                //将PixBlock的值组成数组  
                Array pSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
                double NoDataValue = 0;

                double.TryParse(rasterProps.NoDataValue.ToString(), out NoDataValue);
                Geoprocessor tGp = new Geoprocessor();

                //负向指标

                double tmpMax = 0;
                double tmpMin = 0;

                ESRI.ArcGIS.DataManagementTools.GetRasterProperties GpGetRasterProp = new ESRI.ArcGIS.DataManagementTools.GetRasterProperties(rasterLayer);
                GpGetRasterProp.property_type = "MAXIMUM";
                tGp.Execute(GpGetRasterProp, null);
                tmpMax = (double)GpGetRasterProp.property;
                GpGetRasterProp.property_type = "MINIMUM";
                tGp.Execute(GpGetRasterProp, null);
                tmpMin = (double)GpGetRasterProp.property;
                double length = tmpMax - tmpMin;
                int times = 0;
                //double a1 = tmpMax;
                //double a2 = (tmpMax > (-tmpMin)) ? tmpMax : (-tmpMin);
                #region 极差化
                for (int k = 0; k < rasterProps.Height; k++)
                {
                    for (int j = 0; j < rasterProps.Width; j++)
                    {
                        if (times++ >= 2000)
                        {
                            Application.DoEvents();
                            statusBarXY.Text = statusText + "影像分辨率为:" + rasterProps.Height + "×" + rasterProps.Width + ";正在处理第" + (k + 1) + "行,第" + (j + 1) + "列";
                            times = 0;
                        }
                        double value = 0;
                        double.TryParse(pSafeArray.GetValue(j, k).ToString(), out value);
                        if (Math.Abs(value) > 10000000000000)
                            continue;
                        if (value >= 0)
                            pSafeArray.SetValue((float)(value / tmpMax), j, k);
                        else
                            pSafeArray.SetValue((float)1, j, k);
                    }
                }
                #endregion




                pPixelBlock.set_SafeArray(0, pSafeArray);
                IRasterEdit rasterEdit = pRasters as IRasterEdit;
                rasterEdit.Write(tlp, pPixelBlock);
                rasterEdit.Refresh();
            }
            MessageBox.Show("标准化完成。");
        }

        private void 预警建筑密度_Click(object sender, EventArgs e)
        {

            for (int i = 0; i < axMapControl1.LayerCount; i++)
            {
                IRasterLayer rasterLayer = axMapControl1.get_Layer(i) as IRasterLayer;
                string statusText = "共有：" + axMapControl1.LayerCount + "个图层，正在处理第" + (i + 1) + "个;";
                if (rasterLayer == null)
                    continue;
                IRaster pRasters = rasterLayer.Raster;
                IRaster2 pRaster2s = pRasters as IRaster2;
                IRasterProps rasterProps = (IRasterProps)pRasters;

                IPnt pBlockSize = new Pnt();
                pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);

                double dX = rasterProps.MeanCellSize().X; //栅格的宽度
                double dY = rasterProps.MeanCellSize().Y; //栅格的高度
                //选取整个范围  
                IPixelBlock pPixelBlock = pRasters.CreatePixelBlock(pBlockSize);
                //左上点坐标  
                IPnt tlp = new Pnt();
                tlp.SetCoords(0, 0);
                //读入栅格  
                IRasterBandCollection pRasterBands = pRasters as IRasterBandCollection;
                IRasterBand pRasterBand = pRasterBands.Item(0);
                IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
                pRawRixels.Read(tlp, pPixelBlock);
                //将PixBlock的值组成数组  
                Array pSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
                double NoDataValue = 0;
                int areaSum = 101 * 101;
                int times = 0;
                //double a1 = tmpMax;
                //double a2 = (tmpMax > (-tmpMin)) ? tmpMax : (-tmpMin);
                #region 极差化
                for (int k = 0; k < rasterProps.Height; k++)
                {
                    for (int j = 0; j < rasterProps.Width; j++)
                    {
                        int AreaCount = 0;
                        if (times++ >= 2000)
                        {
                            Application.DoEvents();
                            statusBarXY.Text = statusText + "影像分辨率为:" + rasterProps.Height + "×" + rasterProps.Width + ";正在处理第" + (k + 1) + "行,第" + (j + 1) + "列";
                            times = 0;
                        }
                        double value = 0;
                        double.TryParse(pSafeArray.GetValue(j, k).ToString(), out value);
                        if (Math.Abs(value) > 10000000000000)
                            continue;

                        for (int m = k - 50; m <= k + 50; m++)
                        {
                            for (int n = j - 50; n <= j + 50; n++)
                            {
                                if (m < 0 || n < 0 || n >= rasterProps.Width || m >= rasterProps.Height)
                                    continue;
                                double tmpValue = 0;
                                double.TryParse(pSafeArray.GetValue(n, m).ToString(), out tmpValue);
                                if (Math.Abs(tmpValue) > 10000000000000)
                                    continue;

                                AreaCount++;
                            }
                        }
                        float result = (float)AreaCount / areaSum;
                        pSafeArray.SetValue(result, j, k);
                    }
                }
                #endregion




                pPixelBlock.set_SafeArray(0, pSafeArray);
                IRasterEdit rasterEdit = pRasters as IRasterEdit;
                rasterEdit.Write(tlp, pPixelBlock);
                rasterEdit.Refresh();
            }
            MessageBox.Show("计算完成。");
        }

        private void 人均绿地面积降速ToolStripMenuItem1_Click(object sender, EventArgs e)
        {


            ILayer pointLayer = ArcEngineFunction.GetLayerByName(axMapControl1, "17points");
            IRaster ldRasterLayer12 = (ArcEngineFunction.GetLayerByName(axMapControl1, "17绿地面积.tif") as IRasterLayer).Raster;//12绿地对应的图层
            IRaster ldRasterLayer09 = (ArcEngineFunction.GetLayerByName(axMapControl1, "12绿地面积.tif") as IRasterLayer).Raster;//09绿地对应的图层
            //ILayer popPointLayer = ArcEngineFunction.GetLayerByName(axMapControl1, "12pop_points");
            IRaster pRJRasters = (ArcEngineFunction.GetLayerByName(axMapControl1, "17人均绿地面积降速.tif") as IRasterLayer).Raster;
            IRaster pPOPRasters = (ArcEngineFunction.GetLayerByName(axMapControl1, "17pop.tif") as IRasterLayer).Raster;//人口对应的点图层
            IRaster pPOPRasters2 = (ArcEngineFunction.GetLayerByName(axMapControl1, "12pop.tif") as IRasterLayer).Raster;//人口对应的点图层
            int cellcounts = (pointLayer as IFeatureLayer).FeatureClass.FeatureCount(null);

            #region 人均绿地面积降速
            IRaster2 pRaster2s = pRJRasters as IRaster2;
            IRasterProps rasterProps = (IRasterProps)pRJRasters;
            IPnt pBlockSize = new Pnt();
            pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);

            double dX = rasterProps.MeanCellSize().X; //栅格的宽度
            double dY = rasterProps.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock = pRJRasters.CreatePixelBlock(pBlockSize);
            //左上点坐标  
            IPnt tlp = new Pnt();
            tlp.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands = pRJRasters as IRasterBandCollection;
            IRasterBand pRasterBand = pRasterBands.Item(0);
            IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
            pRawRixels.Read(tlp, pPixelBlock);
            //将PixBlock的值组成数组  
            Array pRJSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
            #endregion
            #region  12年绿地面积
            IRaster2 pRaster2s12 = ldRasterLayer12 as IRaster2;
            IRasterProps rasterProps12 = (IRasterProps)ldRasterLayer12;
            IPnt pBlockSize12 = new Pnt();
            pBlockSize12.SetCoords(rasterProps12.Width, rasterProps12.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock12 = ldRasterLayer12.CreatePixelBlock(pBlockSize12);
            //左上点坐标  
            IPnt tlp12 = new Pnt();
            tlp12.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands12 = ldRasterLayer12 as IRasterBandCollection;
            IRasterBand pRasterBand12 = pRasterBands12.Item(0);
            IRawPixels pRawRixels12 = pRasterBands12.Item(0) as IRawPixels;
            pRawRixels12.Read(tlp12, pPixelBlock12);
            //将PixBlock的值组成数组  
            Array pldSafeArray12 = pPixelBlock12.get_SafeArray(0) as System.Array;
            int LDHeight12 = pldSafeArray12.GetLength(1);
            int LDWidth12 = pldSafeArray12.GetLength(0);
            #endregion

            #region  09年绿地面积
            IRaster2 pRaster2s09 = ldRasterLayer09 as IRaster2;
            IRasterProps rasterProps09 = (IRasterProps)ldRasterLayer09;
            IPnt pBlockSize09 = new Pnt();
            pBlockSize09.SetCoords(rasterProps09.Width, rasterProps09.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock09 = ldRasterLayer09.CreatePixelBlock(pBlockSize09);
            //左上点坐标  
            IPnt tlp09 = new Pnt();
            tlp09.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands09 = ldRasterLayer09 as IRasterBandCollection;
            IRasterBand pRasterBand09 = pRasterBands09.Item(0);
            IRawPixels pRawRixels09 = pRasterBands09.Item(0) as IRawPixels;
            pRawRixels09.Read(tlp09, pPixelBlock09);
            //将PixBlock的值组成数组  
            Array pldSafeArray09 = pPixelBlock09.get_SafeArray(0) as System.Array;
            int LDHeight09 = pldSafeArray09.GetLength(1);
            int LDWidth09 = pldSafeArray09.GetLength(0);
            #endregion
            #region 12人口图层
            IRaster2 pRaster2s1 = pPOPRasters as IRaster2;
            IRasterProps rasterProps1 = (IRasterProps)pPOPRasters;
            IPnt pBlockSize1 = new Pnt();
            pBlockSize1.SetCoords(rasterProps1.Width, rasterProps1.Height);


            double dX1 = rasterProps1.MeanCellSize().X; //栅格的宽度
            double dY1 = rasterProps1.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock1 = pPOPRasters.CreatePixelBlock(pBlockSize1);
            //左上点坐标  
            IPnt tlp1 = new Pnt();
            tlp1.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands1 = pPOPRasters as IRasterBandCollection;
            IRasterBand pRasterBand1 = pRasterBands1.Item(0);
            IRawPixels pRawRixels1 = pRasterBands1.Item(0) as IRawPixels;
            pRawRixels1.Read(tlp1, pPixelBlock1);
            //将PixBlock的值组成数组  
            Array pPOPSafeArray1 = pPixelBlock1.get_SafeArray(0) as System.Array;
            int popHeight = pPOPSafeArray1.GetLength(1);
            int popWidth = pPOPSafeArray1.GetLength(0);
            #endregion

            #region 09人口图层
            IRaster2 pRaster2s2 = pPOPRasters2 as IRaster2;
            IRasterProps rasterProps2 = (IRasterProps)pPOPRasters2;
            IPnt pBlockSize2 = new Pnt();
            pBlockSize2.SetCoords(rasterProps2.Width, rasterProps2.Height);


            double dX2 = rasterProps2.MeanCellSize().X; //栅格的宽度
            double dY2 = rasterProps2.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock2 = pPOPRasters2.CreatePixelBlock(pBlockSize2);
            //左上点坐标  
            IPnt tlp2 = new Pnt();
            tlp2.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands2 = pPOPRasters2 as IRasterBandCollection;
            IRasterBand pRasterBand2 = pRasterBands2.Item(0);
            IRawPixels pRawRixels2 = pRasterBands2.Item(0) as IRawPixels;
            pRawRixels2.Read(tlp2, pPixelBlock2);
            //将PixBlock的值组成数组  
            Array pPOPSafeArray2 = pPixelBlock2.get_SafeArray(0) as System.Array;
            int popHeight2 = pPOPSafeArray2.GetLength(1);
            int popWidth2 = pPOPSafeArray2.GetLength(0);
            #endregion
            string statusBarStr = "";
            int times = 0;
            for (int i = 0; i < cellcounts; i++)
            {
                times++;
                if (times >= 200)
                {
                    Application.DoEvents();
                    statusBarStr = i + "//" + cellcounts;
                    this.statusBarXY.Text = statusBarStr;
                    times = 0;
                }
                IFeature f = (pointLayer as IFeatureLayer).FeatureClass.GetFeature(i);
                IPoint p = f.Shape as IPoint;

                double popSum = 0;
                double popSum09 = 0;
                double LDArea12 = 0;
                double LDArea09 = 0;

                int tmpy = pRaster2s12.ToPixelRow(p.Y);
                int tmpx = pRaster2s12.ToPixelColumn(p.X);
                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= LDWidth12 || k >= LDHeight12)
                            continue;
                        double value = 0;
                        object oo = pldSafeArray12.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value == 0.0 || value == 255)
                            continue;
                        LDArea12++;
                    }
                }
                //LDArea *= 900;

                tmpy = pRaster2s09.ToPixelRow(p.Y);
                tmpx = pRaster2s09.ToPixelColumn(p.X);
                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= LDWidth09 || k >= LDHeight09)
                            continue;
                        double value = 0;
                        object oo = pldSafeArray09.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value == 0.0 || value == 255)
                            continue;
                        LDArea09++;
                    }
                }
                //LDArea *= 900;


                tmpy = pRaster2s1.ToPixelRow(p.Y);
                tmpx = pRaster2s1.ToPixelColumn(p.X);

                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= popWidth || k >= popHeight)
                            continue;
                        double value = 0;
                        object oo = pPOPSafeArray1.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value < -200)
                            continue;
                        popSum += value;
                    }
                }
                popSum *= 0.0009;
                if (popSum < 1)
                    popSum = 1;

                int y = pRaster2s.ToPixelRow(p.Y);
                int x = pRaster2s.ToPixelColumn(p.X);
                if (LDArea09 == 0)
                {
                    float tmp = (float)((-LDArea12 * 900) / popSum);
                    pRJSafeArray.SetValue(tmp, x, y);
                    continue;
                }
                tmpy = pRaster2s2.ToPixelRow(p.Y);
                tmpx = pRaster2s2.ToPixelColumn(p.X);

                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= popWidth2 || k >= popHeight2)
                            continue;
                        double value = 0;
                        object oo = pPOPSafeArray2.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value < -200)
                            continue;
                        popSum09 += value;
                    }
                }
                popSum09 *= 0.0009;
                if (popSum09 < 1)
                    popSum09 = 1;

                double result = (LDArea09 / popSum09 - LDArea12 / popSum) / (LDArea09 / popSum09);

                // double result = buildArea / landCount;
                pRJSafeArray.SetValue((float)result, x, y);
            }
            pPixelBlock.set_SafeArray(0, pRJSafeArray);
            IRasterEdit rasterEdit = pRJRasters as IRasterEdit;
            rasterEdit.Write(tlp, pPixelBlock);
            rasterEdit.Refresh();
            MessageBox.Show("处理完成。");
        }

        private void 人均绿地面积降速ToolStripMenuItem2_Click(object sender, EventArgs e)
        {

            ILayer pointLayer = ArcEngineFunction.GetLayerByName(axMapControl1, "20points");
            IRaster ldRasterLayer12 = (ArcEngineFunction.GetLayerByName(axMapControl1, "20绿地面积.tif") as IRasterLayer).Raster;//12绿地对应的图层
            IRaster ldRasterLayer09 = (ArcEngineFunction.GetLayerByName(axMapControl1, "12绿地面积.tif") as IRasterLayer).Raster;//09绿地对应的图层
            //ILayer popPointLayer = ArcEngineFunction.GetLayerByName(axMapControl1, "12pop_points");
            IRaster pRJRasters = (ArcEngineFunction.GetLayerByName(axMapControl1, "20人均绿地面积降速.tif") as IRasterLayer).Raster;
            IRaster pPOPRasters = (ArcEngineFunction.GetLayerByName(axMapControl1, "20pop.tif") as IRasterLayer).Raster;//人口对应的点图层
            IRaster pPOPRasters2 = (ArcEngineFunction.GetLayerByName(axMapControl1, "12pop.tif") as IRasterLayer).Raster;//人口对应的点图层
            int cellcounts = (pointLayer as IFeatureLayer).FeatureClass.FeatureCount(null);

            #region 人均绿地面积降速
            IRaster2 pRaster2s = pRJRasters as IRaster2;
            IRasterProps rasterProps = (IRasterProps)pRJRasters;
            IPnt pBlockSize = new Pnt();
            pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);

            double dX = rasterProps.MeanCellSize().X; //栅格的宽度
            double dY = rasterProps.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock = pRJRasters.CreatePixelBlock(pBlockSize);
            //左上点坐标  
            IPnt tlp = new Pnt();
            tlp.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands = pRJRasters as IRasterBandCollection;
            IRasterBand pRasterBand = pRasterBands.Item(0);
            IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
            pRawRixels.Read(tlp, pPixelBlock);
            //将PixBlock的值组成数组  
            Array pRJSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
            #endregion
            #region  12年绿地面积
            IRaster2 pRaster2s12 = ldRasterLayer12 as IRaster2;
            IRasterProps rasterProps12 = (IRasterProps)ldRasterLayer12;
            IPnt pBlockSize12 = new Pnt();
            pBlockSize12.SetCoords(rasterProps12.Width, rasterProps12.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock12 = ldRasterLayer12.CreatePixelBlock(pBlockSize12);
            //左上点坐标  
            IPnt tlp12 = new Pnt();
            tlp12.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands12 = ldRasterLayer12 as IRasterBandCollection;
            IRasterBand pRasterBand12 = pRasterBands12.Item(0);
            IRawPixels pRawRixels12 = pRasterBands12.Item(0) as IRawPixels;
            pRawRixels12.Read(tlp12, pPixelBlock12);
            //将PixBlock的值组成数组  
            Array pldSafeArray12 = pPixelBlock12.get_SafeArray(0) as System.Array;
            int LDHeight12 = pldSafeArray12.GetLength(1);
            int LDWidth12 = pldSafeArray12.GetLength(0);
            #endregion

            #region  09年绿地面积
            IRaster2 pRaster2s09 = ldRasterLayer09 as IRaster2;
            IRasterProps rasterProps09 = (IRasterProps)ldRasterLayer09;
            IPnt pBlockSize09 = new Pnt();
            pBlockSize09.SetCoords(rasterProps09.Width, rasterProps09.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock09 = ldRasterLayer09.CreatePixelBlock(pBlockSize09);
            //左上点坐标  
            IPnt tlp09 = new Pnt();
            tlp09.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands09 = ldRasterLayer09 as IRasterBandCollection;
            IRasterBand pRasterBand09 = pRasterBands09.Item(0);
            IRawPixels pRawRixels09 = pRasterBands09.Item(0) as IRawPixels;
            pRawRixels09.Read(tlp09, pPixelBlock09);
            //将PixBlock的值组成数组  
            Array pldSafeArray09 = pPixelBlock09.get_SafeArray(0) as System.Array;
            int LDHeight09 = pldSafeArray09.GetLength(1);
            int LDWidth09 = pldSafeArray09.GetLength(0);
            #endregion
            #region 12人口图层
            IRaster2 pRaster2s1 = pPOPRasters as IRaster2;
            IRasterProps rasterProps1 = (IRasterProps)pPOPRasters;
            IPnt pBlockSize1 = new Pnt();
            pBlockSize1.SetCoords(rasterProps1.Width, rasterProps1.Height);


            double dX1 = rasterProps1.MeanCellSize().X; //栅格的宽度
            double dY1 = rasterProps1.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock1 = pPOPRasters.CreatePixelBlock(pBlockSize1);
            //左上点坐标  
            IPnt tlp1 = new Pnt();
            tlp1.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands1 = pPOPRasters as IRasterBandCollection;
            IRasterBand pRasterBand1 = pRasterBands1.Item(0);
            IRawPixels pRawRixels1 = pRasterBands1.Item(0) as IRawPixels;
            pRawRixels1.Read(tlp1, pPixelBlock1);
            //将PixBlock的值组成数组  
            Array pPOPSafeArray1 = pPixelBlock1.get_SafeArray(0) as System.Array;
            int popHeight = pPOPSafeArray1.GetLength(1);
            int popWidth = pPOPSafeArray1.GetLength(0);
            #endregion

            #region 09人口图层
            IRaster2 pRaster2s2 = pPOPRasters2 as IRaster2;
            IRasterProps rasterProps2 = (IRasterProps)pPOPRasters2;
            IPnt pBlockSize2 = new Pnt();
            pBlockSize2.SetCoords(rasterProps2.Width, rasterProps2.Height);


            double dX2 = rasterProps2.MeanCellSize().X; //栅格的宽度
            double dY2 = rasterProps2.MeanCellSize().Y; //栅格的高度
            //选取整个范围  
            IPixelBlock pPixelBlock2 = pPOPRasters2.CreatePixelBlock(pBlockSize2);
            //左上点坐标  
            IPnt tlp2 = new Pnt();
            tlp2.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands2 = pPOPRasters2 as IRasterBandCollection;
            IRasterBand pRasterBand2 = pRasterBands2.Item(0);
            IRawPixels pRawRixels2 = pRasterBands2.Item(0) as IRawPixels;
            pRawRixels2.Read(tlp2, pPixelBlock2);
            //将PixBlock的值组成数组  
            Array pPOPSafeArray2 = pPixelBlock2.get_SafeArray(0) as System.Array;
            int popHeight2 = pPOPSafeArray2.GetLength(1);
            int popWidth2 = pPOPSafeArray2.GetLength(0);
            #endregion
            string statusBarStr = "";
            int times = 0;
            for (int i = 0; i < cellcounts; i++)
            {
                times++;
                if (times >= 200)
                {
                    Application.DoEvents();
                    statusBarStr = i + "//" + cellcounts;
                    this.statusBarXY.Text = statusBarStr;
                    times = 0;
                }
                IFeature f = (pointLayer as IFeatureLayer).FeatureClass.GetFeature(i);
                IPoint p = f.Shape as IPoint;

                double popSum = 0;
                double popSum09 = 0;
                double LDArea12 = 0;
                double LDArea09 = 0;

                int tmpy = pRaster2s12.ToPixelRow(p.Y);
                int tmpx = pRaster2s12.ToPixelColumn(p.X);
                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= LDWidth12 || k >= LDHeight12)
                            continue;
                        double value = 0;
                        object oo = pldSafeArray12.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value == 0.0 || value == 255)
                            continue;
                        LDArea12++;
                    }
                }
                //LDArea *= 900;

                tmpy = pRaster2s09.ToPixelRow(p.Y);
                tmpx = pRaster2s09.ToPixelColumn(p.X);
                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= LDWidth09 || k >= LDHeight09)
                            continue;
                        double value = 0;
                        object oo = pldSafeArray09.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value == 0.0 || value == 255)
                            continue;
                        LDArea09++;
                    }
                }
                //LDArea *= 900;


                tmpy = pRaster2s1.ToPixelRow(p.Y);
                tmpx = pRaster2s1.ToPixelColumn(p.X);

                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= popWidth || k >= popHeight)
                            continue;
                        double value = 0;
                        object oo = pPOPSafeArray1.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value < -200)
                            continue;
                        popSum += value;
                    }
                }
                popSum *= 0.0009;
                if (popSum < 1)
                    popSum = 1;

                int y = pRaster2s.ToPixelRow(p.Y);
                int x = pRaster2s.ToPixelColumn(p.X);
                if (LDArea09 == 0)
                {
                    float tmp = (float)((-LDArea12 * 900) / popSum);
                    pRJSafeArray.SetValue(tmp, x, y);
                    continue;
                }
                tmpy = pRaster2s2.ToPixelRow(p.Y);
                tmpx = pRaster2s2.ToPixelColumn(p.X);

                for (int k = tmpy - 50; k <= tmpy + 50; k++)
                {
                    for (int j = tmpx - 50; j <= tmpx + 50; j++)
                    {
                        if (k < 0 || j < 0 || j >= popWidth2 || k >= popHeight2)
                            continue;
                        double value = 0;
                        object oo = pPOPSafeArray2.GetValue(j, k);
                        double.TryParse(oo.ToString(), out value);
                        if (value < -200)
                            continue;
                        popSum09 += value;
                    }
                }
                popSum09 *= 0.0009;
                if (popSum09 < 1)
                    popSum09 = 1;

                double result = (LDArea09 / popSum09 - LDArea12 / popSum) / (LDArea09 / popSum09);

                // double result = buildArea / landCount;
                pRJSafeArray.SetValue((float)result, x, y);
            }
            pPixelBlock.set_SafeArray(0, pRJSafeArray);
            IRasterEdit rasterEdit = pRJRasters as IRasterEdit;
            rasterEdit.Write(tlp, pPixelBlock);
            rasterEdit.Refresh();
            MessageBox.Show("处理完成。");
        }

        private void 数据重采样ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            //修复数据不匹配的问题,以第一个栅格图层为准，有空值则赋值为平均值
            //FolderBrowserDialog fbd = new FolderBrowserDialog();
            //if (fbd.ShowDialog() != DialogResult.OK)
            //    return;
            //string savePath = fbd.SelectedPath;


            #region 输出参数
            //int width = rasterProps.Width;
            //int height = rasterProps.Height;
            //double cellsize = rasterProps.MeanCellSize().X;
            //ISpatialReference sr = rasterProps.SpatialReference;
            //IPoint originPoint = rasterProps.Extent.LowerLeft;
            #endregion

            IFeatureLayer pointLayer = axMapControl1.get_Layer(0) as IFeatureLayer;
            int cellcounts = pointLayer.FeatureClass.FeatureCount(null);

            Geoprocessor tGp = new Geoprocessor();
            for (int i = 1; i < axMapControl1.LayerCount; i += 2)
            {
                #region 第一个栅格图层
                IRaster pRasters = (axMapControl1.get_Layer(i + 1) as IRasterLayer).Raster;
                IRaster2 pRaster2s = pRasters as IRaster2;
                IRasterProps rasterProps = (IRasterProps)pRasters;

                IPnt pBlockSize = new Pnt();
                pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);
                //选取整个范围  
                IPixelBlock pPixelBlock = pRasters.CreatePixelBlock(pBlockSize);
                //左上点坐标  
                IPnt tlp = new Pnt();
                tlp.SetCoords(0, 0);
                //读入栅格  
                IRasterBandCollection pRasterBands = pRasters as IRasterBandCollection;
                IRasterBand pRasterBand = pRasterBands.Item(0);
                IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
                pRawRixels.Read(tlp, pPixelBlock);
                //将PixBlock的值组成数组  
                Array pSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
                //object ndv = rasterProps[i - 1].NoDataValue;
                #endregion

                //Array pResult = (Array)pSafeArray.Clone();
                double avg = 0;
                IRasterLayer rasterLayer = axMapControl1.get_Layer(i) as IRasterLayer;
                string statusText = "共有：" + axMapControl1.LayerCount + "个图层，正在处理第" + (i + 1) + "个;";
                if (rasterLayer == null)
                    continue;
                //Normalization(rasterLayer, 1, statusText);

                ESRI.ArcGIS.DataManagementTools.GetRasterProperties GpGetRasterProp = new ESRI.ArcGIS.DataManagementTools.GetRasterProperties(rasterLayer);
                GpGetRasterProp.property_type = "MEAN";
                tGp.Execute(GpGetRasterProp, null);
                avg = (double)GpGetRasterProp.property;

                #region 当前栅格图层
                IRaster pTmpRasters = rasterLayer.Raster;
                IRaster2 pTmpRaster2s = pTmpRasters as IRaster2;
                IRasterProps tmpRasterProps = (IRasterProps)pTmpRasters;

                IPnt pTmpBlockSize = new Pnt();
                pTmpBlockSize.SetCoords(tmpRasterProps.Width, tmpRasterProps.Height);
                //选取整个范围  
                IPixelBlock pTmpPixelBlock = pTmpRasters.CreatePixelBlock(pTmpBlockSize);
                //左上点坐标  
                IPnt tlp1 = new Pnt();
                tlp1.SetCoords(0, 0);
                //读入栅格  
                IRasterBandCollection pTmpRasterBands = pTmpRasters as IRasterBandCollection;
                IRasterBand pTmpRasterBand = pTmpRasterBands.Item(0);
                IRawPixels pTmpRawRixels = pTmpRasterBands.Item(0) as IRawPixels;
                pTmpRawRixels.Read(tlp1, pTmpPixelBlock);
                //将PixBlock的值组成数组  
                Array pTmpSafeArray = pTmpPixelBlock.get_SafeArray(0) as System.Array;
                double NoDataValue = 0;
                double.TryParse(pTmpSafeArray.GetValue(0, 0).ToString(), out NoDataValue);
                #endregion

                int popWidth = pTmpSafeArray.GetLength(0);
                int popHeight = pTmpSafeArray.GetLength(1);
                int times = 0;
                string statusBarStr = "";
                for (int j = 0; j < cellcounts; j++)
                {
                    times++;
                    if (times >= 2000)
                    {
                        Application.DoEvents();
                        statusBarStr = statusText + "   " + j + "/" + cellcounts;
                        this.statusBarXY.Text = statusBarStr;
                        times = 0;
                    }
                    IFeature f = (pointLayer as IFeatureLayer).FeatureClass.GetFeature(j);
                    IPoint p = f.Shape as IPoint;

                    int y = pRaster2s.ToPixelRow(p.Y);
                    int x = pRaster2s.ToPixelColumn(p.X);
                    int tmpy = pTmpRaster2s.ToPixelRow(p.Y);
                    int tmpx = pTmpRaster2s.ToPixelColumn(p.X);
                    if (tmpy < 0 || tmpx < 0 || tmpx >= popWidth || tmpy >= popHeight)
                    {

                        pSafeArray.SetValue((float)avg, x, y);
                        continue;
                    }
                    double value = 0;
                    object oo = pTmpSafeArray.GetValue(tmpx, tmpy);
                    double.TryParse(oo.ToString(), out value);
                    //pSafeArray.SetValue((float)0, x, y);
                    if (value == NoDataValue)
                    {
                        pSafeArray.SetValue((float)avg, x, y);
                        continue;
                    }

                    if (value <= 0)
                    {
                        pSafeArray.SetValue(0, x, y);
                        continue;
                    }

                    pSafeArray.SetValue((float)value, x, y);

                }
                pPixelBlock.set_SafeArray(0, pSafeArray);
                IRasterEdit rasterEdit = pRasters as IRasterEdit;
                rasterEdit.Write(tlp, pPixelBlock);
                rasterEdit.Refresh();

                //statusBarXY.Text = "正在重建信息表。";
                //string LayerPath = (axMapControl1.get_Layer(i + 1) as IRasterLayer).FilePath;
                //statusBarXY.Text = statusText + "正在构建信息表。";
                //ESRI.ArcGIS.DataManagementTools.CalculateStatistics rasterStatistics = new ESRI.ArcGIS.DataManagementTools.CalculateStatistics(LayerPath);
                //tGp.Execute(rasterStatistics, null);
                //IRasterWorkspace2 rws = AEFunction.ArcEngineFunction.CreateRasterWorkspace(savePath);
                //ArcEngineFunction.CreateFileRasterDataset(rws, rasterLayer.Name, originPoint, width, height, cellsize, rstPixelType.PT_FLOAT, sr,pResult)
                //System.Runtime.InteropServices.Marshal.ReleaseComObject(rws);
            }
            MessageBox.Show("数据修复完成！");
        }

        private void 读取txt并进行预警ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            FolderBrowserDialog fbd = new FolderBrowserDialog();
            if (fbd.ShowDialog() != DialogResult.OK)
                return;
            string txtPath = fbd.SelectedPath;

            if (fbd.ShowDialog() != DialogResult.OK)
                return;
            string savePath = fbd.SelectedPath;
            List<string> AllFilePath = GetAllFilePath(txtPath);

            int indicatorCount = AllFilePath.Count;
            ILayer pointLayer = axMapControl1.get_Layer(0);
            int cellcounts = (pointLayer as IFeatureLayer).FeatureClass.FeatureCount(null);
            #region 一个栅格图层
            ILayer layer = axMapControl1.get_Layer(1);
            IRaster pRasters = (layer as IRasterLayer).Raster;
            IRaster2 pRaster2s = pRasters as IRaster2;
            IRasterProps rasterProps = (IRasterProps)pRasters;

            IPnt pBlockSize = new Pnt();
            pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);
            //选取整个范围  
            IPixelBlock pPixelBlock = pRasters.CreatePixelBlock(pBlockSize);
            //左上点坐标  
            IPnt tlp = new Pnt();
            tlp.SetCoords(0, 0);
            //读入栅格  
            IRasterBandCollection pRasterBands = pRasters as IRasterBandCollection;
            IRasterBand pRasterBand = pRasterBands.Item(0);
            IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
            pRawRixels.Read(tlp, pPixelBlock);
            //将PixBlock的值组成数组  
            Array pSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
            #endregion

            #region 输出参数
            int width = rasterProps.Width;
            int height = rasterProps.Height;
            //分项预警指标


            double cellsize = rasterProps.MeanCellSize().X;

            ISpatialReference sr = rasterProps.SpatialReference;
            IPoint originPoint = rasterProps.Extent.LowerLeft;
            #endregion

            Array[] pResult = new Array[indicatorCount];
            Array pTotal = pSafeArray.Clone() as Array;
            List<RoadObject> AllRoad = new List<RoadObject>();//记录每个栅格
            List<WarnLevel> warnLevels = new List<WarnLevel>();//记录每个指标的警限，顺序按EXCEL表来确定
            List<double> values = new List<double>();
            List<double> TotalValues = new List<double>();
            double[] tmpMax = new double[indicatorCount];
            double[] tmpMin = new double[indicatorCount];
            double[] mids = new double[indicatorCount];//记录每个指标的2/3位置的值，顺序按EXCEL表来确定

            List<double> weight = new List<double>();//记录每个指标的权重

            OpenFileDialog weightDlg = new OpenFileDialog();
            weightDlg.Filter = "指标权重文件(*.txt)|*.txt|所有文件(*.*)|*.*";
            if (weightDlg.ShowDialog() == DialogResult.OK)
            {
                StreamReader sr1 = new StreamReader(weightDlg.FileName);
                string line;
                while ((line = sr1.ReadLine()) != null)
                {
                    weight.Add(Convert.ToDouble(line.Split('_')[1]));
                }
            }


            for (int i = 0; i < indicatorCount; i++)
            {
                pResult[i] = pSafeArray.Clone() as Array;
            }
            int count = 0;
            List<string> XY = new List<string>();

            foreach (string txt in AllFilePath)//从txt中读取像元值
            {
                count++;
                StreamReader sw = new StreamReader(txt);
                string stastuStr = "共" + AllFilePath.Count + "个图层，正在读取第" + count + "个。  ";

                int times = 0;
                for (int j = 0; j < cellcounts; j++)
                {
                    times++;
                    if (times >= 2000)
                    {
                        System.Windows.Forms.Application.DoEvents();
                        statusBarXY.Text = stastuStr + (j + 1) + "/" + cellcounts;
                        times = 0;
                    }

                    int y = 0;
                    int x = 0;
                    if (count == 1)
                    {
                        IFeature f = (pointLayer as IFeatureLayer).FeatureClass.GetFeature(j);
                        IPoint p = f.Shape as IPoint;
                        x = pRaster2s.ToPixelColumn(p.X);
                        y = pRaster2s.ToPixelRow(p.Y);
                        XY.Add(x + "_" + y);
                    }
                    else
                    {
                        string[] tmp = XY[j].Split('_');
                        x = Convert.ToInt32(tmp[0]);
                        y = Convert.ToInt32(tmp[1]);
                    }
                    double value = Convert.ToDouble(sw.ReadLine());

                    if (count == 1)
                    {
                        RoadObject ro = new RoadObject(y + "_" + x);
                        ro.Indicators.Add(value);
                        AllRoad.Add(ro);
                    }
                    else
                    {
                        AllRoad[j].Indicators.Add(value);
                    }
                }
                sw.Close();
                GC.Collect();
            }

            int roadCount = AllRoad.Count;

            for (int i = 0; i < AllRoad[0].Indicators.Count; i++)//警限的确定
            {
                foreach (RoadObject ro in AllRoad)
                {
                    values.Add(ro.Indicators[i]);
                }
                values.Sort();
                tmpMin[i] = values[0];
                int position = (int)(roadCount / 3);
                mids[i] = values[position];
                tmpMax[i] = values[roadCount - 1];
                //deltas.Add(stdev(aa.ToArray()));
                values.Clear();
                double[] tmpWL = GetWarnLevelValue(tmpMin[i], mids[i], tmpMax[i]);
                warnLevels.Add(new WarnLevel(tmpWL[0], tmpWL[1], tmpWL[2], tmpWL[3], tmpWL[4], tmpWL[5]));
            }

            int times1 = 0;
            int count1 = 0;
            foreach (RoadObject ro in AllRoad)//计算指标
            {
                count1++;
                times1++;
                if (times1 >= 2000)
                {
                    System.Windows.Forms.Application.DoEvents();
                    statusBarXY.Text = "正在计算第" + count1 + "个指标，共" + cellcounts + "个";
                    times1 = 0;
                }

                string[] X_Y = ro.Name.Split('_');
                int row = Convert.ToInt32(X_Y[0]);
                int col = Convert.ToInt32(X_Y[1]);
                List<double> IndicatorScores = new List<double>();//关联度矩阵
                for (int i = 0; i < indicatorCount; i++)
                {
                    double score = 0; ;

                    score = GetScore(ro.Indicators[i], warnLevels[i]);
                    IndicatorScores.Add(score);

                    int result = GetMaxIndex(score);
                    pResult[i].SetValue((float)result, col, row);
                }
                double totleScore = 0;//综合关联度
                for (int i = 0; i < indicatorCount; i++)
                {
                    totleScore += weight[i] * IndicatorScores[i];
                }
                int totalresult = GetMaxIndex(totleScore);
                pTotal.SetValue((float)totalresult, col, row);
            }

            IRasterWorkspace2 rws = AEFunction.ArcEngineFunction.CreateRasterWorkspace(savePath);
            for (int i = 0; i < indicatorCount; i++)//输出
            {
                GC.Collect();
                string txt = AllFilePath[i];

                string LayerName = txt.Substring(txt.LastIndexOf("\\") + 3, txt.LastIndexOf('.') - txt.LastIndexOf("\\"));
                if (!LayerName.Contains(".tif"))
                {
                    LayerName += ".tif";
                }
                ArcEngineFunction.CreateFileRasterDataset(rws, "预警结果" + LayerName, originPoint, width, height, cellsize, rstPixelType.PT_FLOAT, sr, pResult[i]);
            }
            ArcEngineFunction.CreateFileRasterDataset(rws, "综合预警指标.tif", originPoint, width, height, cellsize, rstPixelType.PT_FLOAT, sr, pTotal);


            System.Runtime.InteropServices.Marshal.ReleaseComObject(rws);
            MessageBox.Show("预警完成。");

        }
        private List<string> GetAllFilePath(string path)
        {
            List<string> AllPath = new List<string>();
            foreach (string directoryName in Directory.GetDirectories(path, "*", SearchOption.TopDirectoryOnly))
            {
                AllPath.AddRange(GetAllFilePath(directoryName));
            }

            foreach (string fileName in Directory.GetFiles(path, "*.txt"))
            {
                AllPath.Add(fileName);
            }
            return AllPath;
        }

        private void ImportDataFromTXT_Click(object sender, EventArgs e)
        {
            OpenFileDialog ofd = new OpenFileDialog();
            ofd.Filter = "TXT(*.txt)|*.txt";
            string txtPath = "";
            string GDBPath = "";
            if (ofd.ShowDialog() != DialogResult.OK)
                return;
            txtPath = ofd.FileName;
            string fileName = ofd.SafeFileName.Split('.')[0];
            FolderBrowserDialog fbd = new FolderBrowserDialog();
            if (fbd.ShowDialog() != DialogResult.OK)
                return;
            GDBPath = fbd.SelectedPath;

            for (int i = 7; i <= 9; i++)
            {
                StreamReader sr = new StreamReader(txtPath, System.Text.Encoding.UTF8);
                string title = sr.ReadLine();
                if (title == "")
                    return;
                string[] fields = title.Split(',');
                fileName += "_" + i;
                CreatePointGDB(ref fileName, GDBPath, fields);
                title = sr.ReadLine();

                IWorkspaceFactory pWSF = new FileGDBWorkspaceFactoryClass();
                IFeatureWorkspace pWS = (IFeatureWorkspace)pWSF.OpenFromFile(GDBPath, 0);
                IFeatureClass newFeatureClass = pWS.OpenFeatureClass(fileName);

                IDataset dataset1 = (IDataset)newFeatureClass;
                IWorkspace workspace1 = dataset1.Workspace;
                //Cast for an IWorkspaceEdit  
                IWorkspaceEdit workspaceEdit1 = (IWorkspaceEdit)workspace1;
                //start an edit session and operation              
                workspaceEdit1.StartEditing(true);
                workspaceEdit1.StartEditOperation();
                int featurecount = 0;
                int times = 0;
                while (title != null)
                {
                    featurecount++;
                    times++;
                    if (times % 5000 == 0)
                    {
                        Application.DoEvents();
                        GC.Collect();
                        statusBarXY.Text = featurecount.ToString();
                        times = 0;
                    }

                    fields = title.Split(',');
                    if (fields[1] != i.ToString())
                    {
                        title = sr.ReadLine();
                        continue;
                    }
                    double x = Convert.ToDouble(fields[fields.Length - 2]);
                    double y = Convert.ToDouble(fields[fields.Length - 1]);
                    if (x < 113.754145 || x > 114.583202 || y < 22.464515 || y > 22.833542)
                    {
                        title = sr.ReadLine();
                        continue;
                    }
                    IPoint p = new PointClass();
                    p.X = x;
                    p.Y = y;
                    try
                    {
                        IFeature f1 = null;
                        f1 = newFeatureClass.CreateFeature();
                        f1.Shape = p;
                        for (int j = 0; j < fields.Length - 2; j++)
                        {
                            f1.set_Value(j + 1, fields[j]);
                        }
                        f1.Store();
                    }
                    catch
                    { }
                    title = sr.ReadLine();
                }

                workspaceEdit1.StopEditOperation();
                workspaceEdit1.StopEditing(true);
                IFeatureClassManage featureClassManage1 = (IFeatureClassManage)newFeatureClass;
                featureClassManage1.UpdateExtent();
                sr.Close();
            }
            MessageBox.Show("导入完成");
        }

        public void CreatePointGDB(ref string strShapeName, string strShapeFolder, string[] fileds)
        {
            //打开工作空间  
            const string strShapeFieldName = "Shape";

            IWorkspaceFactory pWSF = new FileGDBWorkspaceFactoryClass();
            IFeatureWorkspace pWS = (IFeatureWorkspace)pWSF.OpenFromFile(strShapeFolder, 0);

            //设置字段集  
            IFields pFields = new FieldsClass();
            IFieldsEdit pFieldsEdit = (IFieldsEdit)pFields;

            //设置字段  
            IField pField = new FieldClass();
            IFieldEdit pFieldEdit = (IFieldEdit)pField;

            ////创建类型为几何类型的字段  
            pFieldEdit.Name_2 = strShapeFieldName;
            pFieldEdit.Type_2 = esriFieldType.esriFieldTypeGeometry;

            //为esriFieldTypeGeometry类型的字段创建几何定义，包括类型和空间参照   
            IGeometryDef pGeoDef = new GeometryDefClass();     //The geometry definition for the field if IsGeometry is TRUE.  
            IGeometryDefEdit pGeoDefEdit = (IGeometryDefEdit)pGeoDef;
            pGeoDefEdit.GeometryType_2 = esriGeometryType.esriGeometryPoint;
            ISpatialReference pSR = new UnknownCoordinateSystemClass();
            pGeoDefEdit.SpatialReference_2 = pSR;
            pSR.SetDomain(0, 180, 0, 180);

            pFieldEdit.GeometryDef_2 = pGeoDef;
            pFieldsEdit.AddField(pField);

            for (int i = 0; i < fileds.Length - 2; i++)
            {
                string f = fileds[i];
                pField = new FieldClass();
                pFieldEdit = (IFieldEdit)pField;
                pFieldEdit.Name_2 = f;
                pFieldEdit.Type_2 = esriFieldType.esriFieldTypeString;
                pFieldEdit.Precision_2 = 7;//数值精度  
                pFieldEdit.Scale_2 = 3;//小数点位数  
                pFieldsEdit.AddField(pField);
            }

            //创建shapefile  
            try
            {
                pWS.CreateFeatureClass(strShapeName, pFields, null, null, esriFeatureType.esriFTSimple, strShapeFieldName, "");
            }
            catch (System.Exception ex)
            {
                MessageBox.Show(ex.Message);
            }

        }

        private void 公交线路ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            Dictionary<int, List<int>> RouteID_BusID = new Dictionary<int, List<int>>();
            StreamReader sr = new StreamReader("E:\\大数据竞赛\\BUS_DATA\\BUS_DATA\\BUS_ROUTE_DIC.txt");
            string txt = sr.ReadLine();
            if (txt == "")
                return;
            txt = sr.ReadLine();
            while (txt != null)
            {
                string[] array = txt.Split(',');
                int routeID = Convert.ToInt32(array[1]);
                int busID = Convert.ToInt32(array[0]);
                if (!RouteID_BusID.ContainsKey(routeID))
                    RouteID_BusID.Add(routeID, new List<int>());
                RouteID_BusID[routeID].Add(busID);
                txt = sr.ReadLine();
            }
            sr.Close();
            GC.Collect();
            IFeatureLayer pFeatureLayer = axMapControl1.get_Layer(0) as IFeatureLayer;
            IFeatureClass pFeatureClass = pFeatureLayer.FeatureClass;
            foreach (KeyValuePair<int, List<int>> kvp in RouteID_BusID)
            {
                foreach (int value in kvp.Value)
                {

                    IQueryFilter pQueryFilter = new QueryFilterClass();
                    pQueryFilter.WhereClause = string.Format("bus_id='{0}' ", value.ToString());

                    //IQueryFilterDefinition pQueryFilterDefinition = (IQueryFilterDefinition)pQueryFilter;
                    //pQueryFilterDefinition.PostfixClause = "ORDER BY time";

                    DateTime t1 = DateTime.Now;
                    IFeatureCursor pFeatureCursor = pFeatureClass.Search(pQueryFilter, false);
                    List<int> times = new List<int>();
                    IFeature pFeature = pFeatureCursor.NextFeature();
                    while (pFeature != null)
                    {
                        string time = pFeature.get_Value(3) as string;
                        times.Add(Convert.ToInt32(time));
                        pFeature = pFeatureCursor.NextFeature();
                    }
                    DateTime t2 = DateTime.Now;
                    MessageBox.Show((t2 - t1).ToString());
                }
            }
        }

        private void toolStripMenuItem4_Click(object sender, EventArgs e)
        {
            frmClipByTiles clipByTiles = new frmClipByTiles(axMapControl1);
            clipByTiles.ShowDialog();
        }
        /// <summary>
        /// 将多个点图层合并到第一个，并以图层名加入到字段
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void 将多个点图层合并成一个ToolStripMenuItem_Click(object sender, EventArgs e)
        {
            IFeatureLayer netLayer = axMapControl1.get_Layer(0) as IFeatureLayer;
            IFeatureClass netClass = netLayer.FeatureClass;
            IDataset dataset1 = (IDataset)netClass;
            IWorkspace workspace1 = dataset1.Workspace;
            //Cast for an IWorkspaceEdit  
            IWorkspaceEdit workspaceEdit1 = (IWorkspaceEdit)workspace1;
            //start an edit session and operation              
            workspaceEdit1.StartEditing(true);
            workspaceEdit1.StartEditOperation();

            for (int i = 1; i < axMapControl1.LayerCount; i++)
            {
                IFeatureLayer tmpLayer = axMapControl1.get_Layer(i) as IFeatureLayer;
                IFeatureCursor cursor = tmpLayer.FeatureClass.Search(null, false);
                IFeature tmpFeature = cursor.NextFeature();
                while (tmpFeature != null)
                {
                    IFeature f = netClass.CreateFeature();
                    f.Shape = tmpFeature.Shape;
                    f.set_Value(3, tmpLayer.Name);
                    f.set_Value(4, i);
                    f.Store();
                    tmpFeature = cursor.NextFeature();
                }
            }

            workspaceEdit1.StopEditOperation();
            workspaceEdit1.StopEditing(true);
            IFeatureClassManage featureClassManage1 = (IFeatureClassManage)netLayer.FeatureClass;
            featureClassManage1.UpdateExtent();
            MessageBox.Show("聚类完成！");
        }

        private void 读取文件并根据字段分别导出shpToolStripMenuItem_Click(object sender, EventArgs e)
        {
            //字段名称集合，按顺序
            List<string> fields = new List<string>();

            OpenFileDialog ofd = new OpenFileDialog();
            ofd.Filter = "文本文件|*.txt|CSV文件|*.csv|所有文件|*.*";
            if (ofd.ShowDialog() != DialogResult.OK)
                return;
            StreamReader sr = new StreamReader(ofd.FileName);
            string str = sr.ReadLine();
            if (str == "" && str != null)
                return;
            //读取字段名
            string[] aryLine = str.Split(',');
            for (int i = 0; i < aryLine.Length; i++)
            {
                fields.Add(aryLine[i]);
            }

            //选择分类字段
            frmSelectField fsf = new frmSelectField(fields);
            if (fsf.ShowDialog() != DialogResult.OK)
                return;
            List<int> fieldsIndexs = new List<int>();//记录分类字段，坐标字段的字段顺序，顺序是：分类字段、X、Y
            fieldsIndexs.Add(fsf.classfyField);
            fieldsIndexs.Add(fsf.counts);
            fieldsIndexs.Add(fsf.x);
            fieldsIndexs.Add(fsf.y);


            //将数据读取到Datatable中
            str = sr.ReadLine();
            //每个POI按着签到时间整理
            Dictionary<string, POIGroupByTime> POIRecords = new Dictionary<string, POIGroupByTime>();
            while (str != "" && str != null)
            {
                POIGroupByTime poi = new POIGroupByTime();
                aryLine = str.Split(',');

                for (int j = 0; j < aryLine.Length; j++)
                {
                    poi.values.Add(aryLine[j]);
                }
                if (POIRecords.ContainsKey(poi.values[0]))
                    POIRecords[poi.values[0]].CountsPerHour[Convert.ToInt32(poi.values[fieldsIndexs[0]])] = Convert.ToInt32(poi.values[fieldsIndexs[1]]);
                else
                {
                    POIRecords.Add(poi.values[0], poi);
                    POIRecords[poi.values[0]].CountsPerHour[Convert.ToInt32(poi.values[fieldsIndexs[0]])] = Convert.ToInt32(poi.values[fieldsIndexs[1]]);
                }

                str = sr.ReadLine();
            }
            sr.Close();


            string fileName = System.IO.Path.GetFileNameWithoutExtension(ofd.FileName);
            string filePath = System.IO.Path.GetDirectoryName(ofd.FileName);
            List<ExcelField> excelField = new List<ExcelField>();
            for (int i = 0; i < fields.Count; i++)
            {
                if (fieldsIndexs.Contains(i))
                    continue;
                ExcelField ef = new ExcelField();
                ef.field = fields[i];
                ef.type = typeof(string);
                excelField.Add(ef);
            }
            for (int i = 0; i < 24; i++)
            {
                ExcelField ef = new ExcelField();
                ef.field = i.ToString();
                ef.type = typeof(int);
                excelField.Add(ef);
            }
            ArcEngineFunction.CreatePointShapefile(ref fileName, filePath, null, excelField);

            axMapControl1.AddShapeFile(filePath, fileName + ".shp");
            ILayer newLayer = axMapControl1.get_Layer(0);
            IFeatureClass fc = AEFunction.ArcEngineFunction.GetLayerFeatureClass(newLayer);
            IDataset dataset = (IDataset)fc;
            IWorkspace workspace = dataset.Workspace;
            //Cast for an IWorkspaceEdit  
            IWorkspaceEdit workspaceEdit = (IWorkspaceEdit)workspace;
            //start an edit session and operation              
            workspaceEdit.StartEditing(true);
            workspaceEdit.StartEditOperation();

            int times = 1;
            string statusTxt = "共有" + POIRecords.Count + "个要素，正在创建第";

            foreach (KeyValuePair<string, POIGroupByTime> kvp in POIRecords)
            {
                Application.DoEvents();
                this.statusBarXY.Text = statusTxt + times + "个";
                POIGroupByTime poi = kvp.Value;
                IFeature f = fc.CreateFeature();
                IPoint p = new PointClass();
                p.X = Convert.ToDouble(poi.values[fieldsIndexs[2]]);
                p.Y = Convert.ToDouble(poi.values[fieldsIndexs[3]]);
                f.Shape = p;
                for (int i = 0, m = 0; i < poi.values.Count; i++)
                {
                    if (fieldsIndexs.Contains(i))
                    {
                        m++;
                        continue;
                    }
                    f.set_Value(i + 2 - m, poi.values[i]);
                }
                int tmp = fields.Count - fieldsIndexs.Count;
                for (int i = 0; i < 24; i++)
                {
                    f.set_Value(tmp + 2 + i, poi.CountsPerHour[i]);
                }
                f.Store();
                times++;
            }
            GC.Collect();
            this.statusBarXY.Text = "创建完成，正在刷新数据。";
            Application.DoEvents();
            workspaceEdit.StopEditOperation();
            workspaceEdit.StopEditing(true);
            IFeatureClassManage featureClassManage = (IFeatureClassManage)fc;
            featureClassManage.UpdateExtent();
            axMapControl1.Extent = axMapControl1.FullExtent;
            MessageBox.Show("要素生成完成。", "提示", MessageBoxButtons.OK);
        }

        private void CalcFeatureCounts_Click(object sender, EventArgs e)
        {
            frmCalcFeatureCounts fcc = new frmCalcFeatureCounts();
            if (fcc.ShowDialog() != DialogResult.OK)
                return;
            int pointFeatureCount = 0;
            int pointVextexCount = 0;
            int lineFeatureCount = 0;
            int lineVextexCount = 0;
            int areaFeatureCount = 0;
            int areaVextexCount = 0;
            string folder = fcc.FileSource;
            string filepath = fcc.resultPaht;
            StreamWriter sw = new StreamWriter(filepath);
            DirectoryInfo dr = new DirectoryInfo(folder);
            FileInfo[] files = dr.GetFiles("*.dbf");
            toolStripProgressBar1.Maximum = files.Length;
            toolStripProgressBar1.Value = 0;
            toolStripProgressBar1.Step = 1;
            foreach (FileInfo f in files)
            {
                toolStripProgressBar1.Value++;
                string tmpName = f.Name.Split('.')[0];
                IWorkspaceFactory pWorkspaceFactory = new ShapefileWorkspaceFactory();


                IWorkspace pWorkspace = pWorkspaceFactory.OpenFromFile(folder, 0);
                IFeatureWorkspace pFeatureWorkspace = pWorkspace as IFeatureWorkspace;
                IFeatureClass pFeatureClass = pFeatureWorkspace.OpenFeatureClass(tmpName + ".shp");

                
                if(pFeatureClass.ShapeType.ToString().Contains("Point"))
                {
                    pointFeatureCount+=pFeatureClass.FeatureCount(null);
                    
                }
                else if((pFeatureClass.ShapeType.ToString().Contains("line")))
                {
                    lineFeatureCount += pFeatureClass.FeatureCount(null);
                }
                else if((pFeatureClass.ShapeType.ToString().Contains("Polygon")))
                {
                    areaFeatureCount += pFeatureClass.FeatureCount(null);
                }
                IFeatureCursor featureCursor = pFeatureClass.Search(null, false);
                IFeature feautre = featureCursor.NextFeature();
                while(feautre!=null)
                {

                    IPointCollection pc = feautre.Shape as IPointCollection;
                    if (pFeatureClass.ShapeType.ToString().Contains("Point"))
                    {

                        pointVextexCount++;

                    }
                    else if ((pFeatureClass.ShapeType.ToString().Contains("line")))
                    {
                        lineVextexCount += pc.PointCount;
                    }
                    else if ((pFeatureClass.ShapeType.ToString().Contains("Polygon")))
                    {
                        areaVextexCount += pc.PointCount;
                    }
                    feautre = featureCursor.NextFeature();
                }
            }
            string str = "{0}要素：{1}个，节点总数：{2}。";
            sw.WriteLine(string.Format(str,"点",pointFeatureCount,pointVextexCount));
            sw.WriteLine(string.Format(str, "线", lineFeatureCount, lineVextexCount));
            sw.WriteLine(string.Format(str, "面", areaFeatureCount, areaVextexCount));
            sw.WriteLine(string.Format(str, "", pointFeatureCount + lineFeatureCount + areaFeatureCount, pointVextexCount + lineVextexCount + areaVextexCount));
            sw.Close();

            toolStripProgressBar1.Maximum = 100;
            toolStripProgressBar1.Value = 0;
            MessageBox.Show("统计完成");
        }



    }

}